time_boxmuller_cuda.cu (6109B)
1 // Test for boxmuller.h on CUDA 2 #include <Random123/philox.h> 3 #include <Random123/threefry.h> 4 #include "util.h" // for timer() 5 #include "util_cuda.h" // for cuda_init, CHECKCALL 6 #include <Random123/boxmuller.hpp> 7 8 typedef r123::Philox4x32 CBRNGF; 9 typedef r123::Threefry2x64 CBRNGD; 10 int debug = 0; 11 const char *progname = "time_boxmuller_cuda"; 12 13 // Sometimes warnings are A LOT more trouble than they're worth. 14 // if we just write u[6], we get warnings 15 // so we write u[(csize>n)?6:0]. 16 #define UGLY(n) (csize>n)?n:0 17 18 // The timedloop kernel sums N randoms per thread for timing and 19 // records that sum in out[tid] (mainly to ensure that 20 // the random generation process does not get optimized away) 21 template <typename CBRNG, typename F, typename F2> 22 __global__ void timedloop(F *out, typename CBRNG::ukey_type k, size_t N){ 23 unsigned tid = blockDim.x * blockIdx.x + threadIdx.x; 24 const size_t klast = sizeof(k)/sizeof(k[0]) - 1; 25 R123_ASSERT(k[klast] == 0); // uses last element of key to 26 k[klast] = tid; // ensure unique key per thread 27 F sum = 0.f; 28 typename CBRNG::ctr_type ctr = {}; 29 const size_t csize = sizeof(ctr)/sizeof(ctr[0]); 30 CBRNG rng; 31 32 for(size_t i=0; i<N; i+=csize){ 33 ctr.incr(); 34 typename CBRNG::ctr_type u = rng(ctr, k); 35 F2 f2; 36 // Using a loop instead of the Duff device here costs 10%, 37 // at least in CUDA4.2 circa Jan 2013 on a Tesla C2050! 38 switch(csize) { 39 case 8: f2 = r123::boxmuller(u[UGLY(6)], u[UGLY(7)]); sum += f2.x + f2.y; 40 f2 = r123::boxmuller(u[UGLY(4)], u[UGLY(5)]); sum += f2.x + f2.y; 41 case 4: f2 = r123::boxmuller(u[UGLY(2)], u[UGLY(3)]); sum += f2.x + f2.y; 42 case 2: f2 = r123::boxmuller(u[0], u[1]); sum += f2.x + f2.y; 43 break; 44 default: 45 R123_ASSERT(0); 46 } 47 } 48 out[tid] = sum; 49 } 50 51 // The dumploop kernel records all the normal randoms individually in out, 52 // so it produces N randoms per thread. Each thread records 53 // its randoms in tid, NTHREADS+tid, NTHREAD*2+tid, ..., NTHREADS*(N-1)+tid 54 // which hopefully results in nicely coalesced writes from each warp. 55 template <typename CBRNG, typename F, typename F2> 56 __global__ void dumploop(F *out, typename CBRNG::ukey_type k, size_t N){ 57 unsigned tid = blockDim.x * blockIdx.x + threadIdx.x; 58 const size_t klast = sizeof(k)/sizeof(k[0]) - 1; 59 R123_ASSERT(k[klast] == 0); // uses last element of key to 60 k[klast] = tid; // ensure unique key per thread 61 typename CBRNG::ctr_type ctr = {}; 62 const size_t csize = sizeof(ctr)/sizeof(ctr[0]); 63 CBRNG rng; 64 65 for(size_t i=0; i<N;){ 66 ctr.incr(); 67 typename CBRNG::ctr_type u = rng(ctr, k); 68 F2 f2; 69 // Using a loop instead of the Duff device here costs 10%, 70 // at least in CUDA4.2 circa Jan 2013 on a Tesla C2050! 71 switch(csize) { 72 case 8: f2 = r123::boxmuller(u[UGLY(6)], u[UGLY(7)]); 73 out[blockDim.x*gridDim.x*i + tid] = f2.x; 74 i++; 75 out[blockDim.x*gridDim.x*i + tid] = f2.y; 76 i++; 77 f2 = r123::boxmuller(u[UGLY(4)], u[UGLY(5)]); 78 out[blockDim.x*gridDim.x*i + tid] = f2.x; 79 i++; 80 out[blockDim.x*gridDim.x*i + tid] = f2.y; 81 i++; 82 case 4: f2 = r123::boxmuller(u[UGLY(2)], u[UGLY(3)]); 83 #undef UGLY 84 out[blockDim.x*gridDim.x*i + tid] = f2.x; 85 i++; 86 out[blockDim.x*gridDim.x*i + tid] = f2.y; 87 i++; 88 case 2: f2 = r123::boxmuller(u[0], u[1]); 89 out[blockDim.x*gridDim.x*i + tid] = f2.x; 90 i++; 91 out[blockDim.x*gridDim.x*i + tid] = f2.y; 92 i++; 93 break; 94 default: 95 asm("trap;"); 96 } 97 } 98 } 99 100 template <typename CBRNG, typename F, typename F2> 101 void timedcall(const char *tname, const char *out_fname, CUDAInfo *infop, typename CBRNG::ukey_type k, size_t N) { 102 double cur_time, dt; 103 const int nthreads = infop->blocks_per_grid*infop->threads_per_block; 104 const size_t nrand = out_fname ? N * nthreads : nthreads; 105 const size_t out_size = nrand*sizeof(F); 106 F *d_out, *h_out = (F *) malloc(out_size); 107 CHECKNOTZERO(h_out); 108 CHECKCALL(cudaMalloc(&d_out, out_size)); 109 (void) timer(&cur_time); 110 if (out_fname) 111 dumploop<CBRNG,F,F2> <<<infop->blocks_per_grid, infop->threads_per_block>>> (d_out, k, N); 112 else 113 timedloop<CBRNG,F,F2> <<<infop->blocks_per_grid, infop->threads_per_block>>> (d_out, k, N); 114 CHECKCALL(cudaDeviceSynchronize()); 115 CHECKCALL(cudaMemcpy(h_out, d_out, out_size, cudaMemcpyDeviceToHost)); 116 dt = timer(&cur_time); 117 printf("%s %zd in %g sec: %gM/sec\n", tname, N*nthreads, dt, N*nthreads*1.e-6/dt); 118 if (out_fname) { 119 char *fname = (char *) malloc(strlen(out_fname) + strlen(tname) + 2); 120 CHECKNOTZERO(fname); 121 sprintf(fname, "%s-%s", out_fname, tname); 122 FILE *fp = fopen(fname, "w"); 123 CHECKNOTZERO(fp); 124 for (size_t i = 0; i < nrand; i++){ 125 fprintf(fp, "%g\n", h_out[i]); 126 } 127 fclose(fp); 128 free(fname); 129 } else { 130 int nwoops = 0; 131 printf("%s h_out[0] = %g\n", tname, h_out[0]); 132 for (size_t i = 0; i < nrand; i++){ 133 if(h_out[i] == 0.f){ 134 if(nwoops++<10) 135 printf("Woops %s h_out[%zd] = %g\n", tname, i, h_out[i]); 136 137 } 138 } 139 if(nwoops>10){ 140 printf("Woops %s %d times\n", tname, nwoops); 141 } 142 } 143 CHECKCALL(cudaFree(d_out)); 144 free(h_out); 145 } 146 147 const size_t DEF_N = 200000; 148 149 int main(int argc, char **argv){ 150 CBRNGF::ukey_type keyf = {}; 151 CBRNGD::ukey_type keyd = {}; 152 size_t Ntry = DEF_N; 153 char *cp = getenv("R123_DEBUG"); 154 if (cp) 155 debug = atoi(cp); 156 if ((cp = getenv("BOXMULLER_DUMPFILE")) != NULL) { 157 Ntry = 8; 158 } else { 159 Ntry = DEF_N; 160 } 161 if(argc>1) { 162 if (argv[1][0] == '-') { 163 fprintf(stderr, "Usage: %s [iterations_per_thread [key0 [key1]]]\n", argv[0]); 164 exit(1); 165 } 166 Ntry = atol(argv[1]); 167 } 168 // XXX cannot use keyf.size in host code, only in device code 169 for (int i = 0; i < (int)(sizeof(keyf)/sizeof(keyf[0])-1) && 2+i < argc; i++) { 170 keyf.v[i] = atol(argv[2+i]); 171 } 172 for (int i = 0; i < (int)(sizeof(keyd)/sizeof(keyd[0])-1) && 2+i < argc; i++) { 173 keyd.v[i] = atol(argv[2+i]); 174 } 175 CUDAInfo *infop = cuda_init(NULL); 176 timedcall<CBRNGF,float,r123::float2>("float", cp, infop, keyf, Ntry); 177 timedcall<CBRNGD,double,r123::double2>("double",cp, infop, keyd, Ntry); 178 cuda_done(infop); 179 return 0; 180 } 181