time_boxmuller.cpp (4137B)
1 // Test for boxmuller.h on CPU 2 #include <Random123/philox.h> 3 #include <Random123/threefry.h> 4 #include <Random123/boxmuller.hpp> 5 #include "util.h" // for timer() 6 7 #if __GNUC__>=7 8 #pragma GCC diagnostic ignored "-Wimplicit-fallthrough" 9 #endif 10 11 typedef r123::Philox4x32 CBRNGF; 12 #if R123_USE_64BIT 13 typedef r123::Threefry2x64 CBRNGD; 14 #endif 15 16 const char *progname = "time_boxmuller"; 17 18 // Each call to boxmuller() returns a pair of values in the .x and .y 19 // members, which we add up into sum just to avoid being optimized away. 20 template <typename CBRNG, typename F, typename F2> 21 F timedloop(typename CBRNG::ukey_type k, size_t Ntry){ 22 F sum = 0.f; 23 typename CBRNG::ctr_type ctr = {{}}; 24 const size_t csize = sizeof(ctr)/sizeof(ctr[0]); 25 CBRNG rng; 26 27 for(size_t i=0; i<Ntry; i+=csize){ 28 ctr.incr(); 29 typename CBRNG::ctr_type u = rng(ctr, k); 30 F2 f2; 31 switch(csize) { 32 case 8: f2 = r123::boxmuller(u[6], u[7]); sum += f2.x + f2.y; 33 f2 = r123::boxmuller(u[4], u[5]); sum += f2.x + f2.y; 34 case 4: f2 = r123::boxmuller(u[2], u[3]); sum += f2.x + f2.y; 35 case 2: f2 = r123::boxmuller(u[0], u[1]); sum += f2.x + f2.y; 36 break; 37 default: 38 R123_ASSERT(0); 39 } 40 } 41 return sum; 42 } 43 44 template <typename CBRNG, typename F2> 45 void dumploop(FILE *fp, typename CBRNG::ukey_type k, size_t Ntry){ 46 typename CBRNG::ctr_type ctr = {{}}; 47 const size_t csize = sizeof(ctr)/sizeof(ctr[0]); 48 CBRNG rng; 49 50 for(size_t i=0; i<Ntry; i+=csize){ 51 ctr.incr(); 52 typename CBRNG::ctr_type u = rng(ctr, k); 53 F2 f2; 54 switch(csize) { 55 case 8: f2 = r123::boxmuller(u[6], u[7]); fprintf(fp, "%g\n%g\n", f2.x, f2.y); 56 f2 = r123::boxmuller(u[4], u[5]); fprintf(fp, "%g\n%g\n", f2.x, f2.y); 57 case 4: f2 = r123::boxmuller(u[2], u[3]); fprintf(fp, "%g\n%g\n", f2.x, f2.y); 58 case 2: f2 = r123::boxmuller(u[0], u[1]); fprintf(fp, "%g\n%g\n", f2.x, f2.y); break; 59 default: 60 R123_ASSERT(0); 61 } 62 } 63 } 64 65 #define NREPEAT 20 66 67 template <typename CBRNG, typename F, typename F2> 68 void timedcall(const char *tname, typename CBRNG::ukey_type k, size_t Ntry, char *out_fname) { 69 double cur_time, dt; 70 F sums[NREPEAT]; 71 int i; 72 FILE *fp; 73 char *fname; 74 if (out_fname) { 75 fname = (char *) malloc(strlen(out_fname) + strlen(tname) + 2); 76 CHECKNOTZERO(fname); 77 sprintf(fname, "%s-%s", out_fname, tname); 78 fp = fopen(fname, "w"); 79 CHECKNOTZERO(fp); 80 } else { 81 fname = NULL; 82 fp = NULL; 83 } 84 (void) timer(&cur_time); 85 /* 86 * we call timedloop NREPEAT times so that it is easy to keep 87 * Ntry the same for boxmuller.cu and boxmuller.cpp, so sum[0] 88 * can be checked. 89 */ 90 for (i = 0; i < NREPEAT; i++) { 91 k.v[sizeof(k)/sizeof(k.v[0])-1] = i; 92 if (fp) 93 dumploop<CBRNG, F2>(fp, k, Ntry); 94 else 95 sums[i] = timedloop<CBRNG, F, F2>(k, Ntry); 96 } 97 dt = timer(&cur_time); 98 if (fp) { 99 printf("%s %lu written to %s in %g sec: %gM/sec\n", tname, (unsigned long)(Ntry*NREPEAT), fname, dt, Ntry*NREPEAT*1.e-6/dt); 100 fclose(fp); 101 free(fname); 102 } else { 103 printf("%s %lu in %g sec: %gM/sec, sum = %g\n", tname, (unsigned long)(Ntry*NREPEAT), dt, Ntry*NREPEAT*1.e-6/dt, sums[0]); 104 for (i = 1; i < NREPEAT; i++) { 105 printf(" %g", sums[i]); 106 } 107 printf("\n"); 108 } 109 } 110 111 const size_t DEF_N = 200000; 112 113 int main(int argc, char **argv){ 114 CBRNGF::ukey_type keyf = {{}}; 115 #if R123_USE_64BIT 116 CBRNGD::ukey_type keyd = {{}}; 117 #endif 118 size_t Ntry = DEF_N; 119 char *dumpfname; 120 121 dumpfname = getenv("BOXMULLER_DUMPFILE"); 122 if(argc>1) { 123 if (argv[1][0] == '-') { 124 fprintf(stderr, "Usage: %s [iterations_per_thread [key0 [key1]]]\n", argv[0]); 125 exit(1); 126 } 127 Ntry = atol(argv[1]); 128 } 129 for (int i = 0; i < (int)(sizeof(keyf)/sizeof(keyf[0])-1) && 2+i < argc; i++) { 130 keyf.v[i] = atol(argv[2+i]); 131 } 132 timedcall<CBRNGF,float,r123::float2>("float", keyf, Ntry, dumpfname); 133 134 #if R123_USE_64BIT 135 for (int i = 0; i < (int)(sizeof(keyd)/sizeof(keyd[0])-1) && 2+i < argc; i++) { 136 keyd.v[i] = atol(argv[2+i]); 137 } 138 timedcall<CBRNGD,double,r123::double2>("double", keyd, Ntry, dumpfname); 139 #endif 140 return 0; 141 } 142