eic-opticks/qudarap/QEvt.cu

0001 #include <stdio.h>
0002
0003 #include "scuda.h"
0004 #include "squad.h"
0005 #include "srec.h"
0006 #include "sphoton.h"
0007 #include "sevent.h"
0008
0009 #include "iexpand.h"
0010 #include "strided_range.h"
0011 #include <thrust/device_vector.h>
0012
0013 /**
0014 _QEvt_checkEvt
0015 -----------------
0016
0017 Demonstrates using seed buffer to lookup genstep_id from photon_id
0018
0019 **/
0020
0021 __global__ void _QEvt_checkEvt(sevent* evt, unsigned width, unsigned height)
0022 {
0023     unsigned ix = blockIdx.x*blockDim.x + threadIdx.x;
0024     if( ix >= width ) return ;
0025
0026     unsigned photon_id = ix ;
0027     unsigned genstep_id = evt->seed[photon_id] ;
0028     const quad6& gs = evt->genstep[genstep_id] ;
0029     int gencode = gs.q0.i.x ;
0030     unsigned num_photon = evt->num_photon ;
0031
0032     printf("//_QEvt_checkEvt width %d height %d photon_id %3d genstep_id %3d  gs.q0.i ( %3d %3d %3d %3d )  gencode %d num_photon %d \n",
0033        width,
0034        height,
0035        photon_id,
0036        genstep_id,
0037        gs.q0.i.x,
0038        gs.q0.i.y,
0039        gs.q0.i.z,
0040        gs.q0.i.w,
0041        gencode,
0042        num_photon
0043       );
0044 }
0045
0046 extern "C" void QEvt_checkEvt(dim3 numBlocks, dim3 threadsPerBlock, sevent* evt, unsigned width, unsigned height )
0047 {
0048     printf("//QEvt_checkEvt width %d height %d \n", width, height );
0049     _QEvt_checkEvt<<<numBlocks,threadsPerBlock>>>( evt, width, height  );
0050 }
0051
0052 /**
0053 QEvt_count_genstep_photons
0054 -------------------------------
0055
0056 NB this needs nvcc compilation due to the use of thrust but
0057 the method itself does not run on the device although the
0058 methods it invokes do run on the device.
0059
0060 So the sevent* argument must be the CPU side instance
0061 which must be is holding GPU side pointers.
0062
0063 **/
0064
0065
0066 //#ifdef DEBUG_QEVENT
0067 struct printf_functor
0068 {
0069     __host__ __device__ void operator()(int x){ printf("printf_functor %d\n", x); }
0070 };
0071 //#endif
0072
0073
0074 /**
0075 QEvt_count_genstep_photons
0076 -----------------------------
0077
0078 Notice how using strided_range needs itemsize stride twice,
0079 because are grabbing single ints "numphoton" from each quad6 6*4 genstep
0080
0081 **/
0082
0083
0084 extern "C" unsigned QEvt_count_genstep_photons(sevent* evt)
0085 {
0086     typedef typename thrust::device_vector<int>::iterator Iterator;
0087
0088     thrust::device_ptr<int> t_gs = thrust::device_pointer_cast( (int*)evt->genstep ) ;
0089
0090 #ifdef DEBUG_QEVENT
0091     printf("//QEvt_count_genstep_photons sevent::genstep_numphoton_offset %d  sevent::genstep_itemsize  %d  \n",
0092             sevent::genstep_numphoton_offset, sevent::genstep_itemsize );
0093 #endif
0094
0095     strided_range<Iterator> gs_pho(
0096         t_gs + sevent::genstep_numphoton_offset,
0097         t_gs + evt->num_genstep*sevent::genstep_itemsize ,
0098         sevent::genstep_itemsize );    // begin, end, stride
0099
0100     evt->num_seed = thrust::reduce(gs_pho.begin(), gs_pho.end() );
0101
0102 #ifdef DEBUG_QEVENT
0103     //thrust::for_each( gs_pho.begin(), gs_pho.end(), printf_functor() );
0104     printf("//QEvt_count_genstep_photons evt.num_genstep %d evt.num_seed %d evt.max_photon %d \n", evt->num_genstep, evt->num_seed, evt->max_photon );
0105 #endif
0106     assert( evt->num_seed <= evt->max_photon );
0107
0108     return evt->num_seed ;
0109 }
0110
0111 /**
0112 QEvt_fill_seed_buffer
0113 -------------------------
0114
0115 Populates seed buffer using the numbers of photons per genstep from the genstep buffer.
0116
0117 See thrustrap/tests/iexpand_stridedTest.cu for the lead up to this
0118
0119 1. use GPU side genstep array to add the numbers of photons
0120    from each genstep giving the total number of photons and seeds *num_seeds*
0121    from all the gensteps
0122
0123 2. populate it by repeating genstep indices into it,
0124    according to the number of photons in each genstep
0125
0126 t_gs+sevent::genstep_numphoton_offset
0127    q0.u.w of the quad6 genstep, which contains the number of photons
0128    for this genstep
0129
0130
0131 WARNING : SOMETHING HERE MESSES UP UNLESS THE SEED BUFFER IS ZEROED PRIOR TO THIS BEING CALLED
0132
0133 ACTUALLY THIS IS DUE TO A A LIMITATION OF IEXPAND, see sysrap/iexpand.h::
0134
0135     NB the output device must be zeroed prior to calling iexpand.
0136     This is because the iexpand is implemented ending with an inclusive_scan
0137     to fill in the non-transition values which relies on initial zeroing.
0138
0139 **/
0140
0141 extern "C" void QEvt_fill_seed_buffer(sevent* evt )
0142 {
0143 #ifdef DEBUG_QEVENT
0144     printf("//QEvt_fill_seed_buffer evt.num_genstep %d evt.num_seed %d evt.max_photon %d \n", evt->num_genstep, evt->num_seed, evt->max_photon );
0145 #endif
0146
0147     assert( evt->seed && evt->num_seed > 0 );
0148     assert( evt->num_seed <= evt->max_photon );
0149
0150     thrust::device_ptr<int> t_seed = thrust::device_pointer_cast(evt->seed) ;
0151
0152     typedef typename thrust::device_vector<int>::iterator Iterator;
0153
0154     thrust::device_ptr<int> t_gs = thrust::device_pointer_cast( (int*)evt->genstep ) ;
0155
0156     strided_range<Iterator> gs_pho(
0157            t_gs + sevent::genstep_numphoton_offset,
0158            t_gs + evt->num_genstep*sevent::genstep_itemsize,
0159            sevent::genstep_itemsize );    // begin, end, stride
0160
0161
0162     //thrust::for_each( gs_pho.begin(), gs_pho.end(), printf_functor() );
0163
0164     iexpand( gs_pho.begin(), gs_pho.end(), t_seed, t_seed + evt->num_seed );
0165
0166     //thrust::for_each( t_seed,  t_seed + evt->num_seed, printf_functor() );
0167
0168 }
0169
0170
0171
0172 /**
0173 QEvt_count_genstep_photons_and_fill_seed_buffer
0174 ---------------------------------------------------
0175
0176 This function does the same as the above two functions.
0177 It is invoked from QEvt::setGenstep
0178
0179 **/
0180
0181 extern "C" void QEvt_count_genstep_photons_and_fill_seed_buffer(sevent* evt )
0182 {
0183     typedef typename thrust::device_vector<int>::iterator Iterator;
0184
0185     thrust::device_ptr<int> t_gs = thrust::device_pointer_cast( (int*)evt->genstep ) ;
0186
0187 #ifdef DEBUG_QEVENT
0188     printf("//QEvt_count_genstep_photons sevent::genstep_numphoton_offset %d  sevent::genstep_itemsize  %d  \n",
0189             sevent::genstep_numphoton_offset, sevent::genstep_itemsize );
0190 #endif
0191
0192
0193     strided_range<Iterator> gs_pho(
0194         t_gs + sevent::genstep_numphoton_offset,
0195         t_gs + evt->num_genstep*sevent::genstep_itemsize ,
0196         sevent::genstep_itemsize );    // begin, end, stride
0197
0198     evt->num_seed = thrust::reduce(gs_pho.begin(), gs_pho.end() );
0199
0200 #ifdef DEBUG_QEVENT
0201     printf("//QEvt_count_genstep_photons_and_fill_seed_buffer evt.num_genstep %d evt.num_seed %d evt.max_photon %d \n", evt->num_genstep, evt->num_seed, evt->max_photon );
0202 #endif
0203
0204     bool expect_seed =  evt->seed && evt->num_seed > 0 ;
0205     if(!expect_seed) printf("//QEvt_count_genstep_photons_and_fill_seed_buffer  evt.seed %s  evt.num_seed %d \n",  (evt->seed ? "YES" : "NO " ), evt->num_seed );
0206     assert( expect_seed );
0207
0208     bool num_seed_ok = evt->num_seed <= evt->max_photon ;
0209
0210     if( num_seed_ok == false )
0211     {
0212         printf("//QEvt_count_genstep_photons_and_fill_seed_buffer FAIL evt.num_seed %d evt.max_photon %d num_seed_ok %d \n", evt->num_seed, evt->max_photon, num_seed_ok  );
0213     }
0214
0215     assert( num_seed_ok );
0216
0217     thrust::device_ptr<int> t_seed = thrust::device_pointer_cast(evt->seed) ;
0218
0219     //thrust::for_each( gs_pho.begin(), gs_pho.end(), printf_functor() );
0220
0221 #ifdef DEBUG_QEVENT
0222     printf("//[QEvt_count_genstep_photons_and_fill_seed_buffer iexpand \n" );
0223 #endif
0224
0225
0226     iexpand( gs_pho.begin(), gs_pho.end(), t_seed, t_seed + evt->num_seed );
0227
0228     //thrust::for_each( t_seed,  t_seed + evt->num_seed, printf_functor() );
0229
0230
0231 #ifdef DEBUG_QEVENT
0232     printf("//]QEvt_count_genstep_photons_and_fill_seed_buffer iexpand \n" );
0233 #endif
0234
0235
0236
0237 }
0238
0239
0240
0241