wyj
11e703e33d
Archive fixed catalog subsets and captured events, manifests, reproducible commands, raw timing and regression logs, and device resource evidence. Document distribution-dependent tile reduction gains and independent indexed producer throughput. Update the GPU plan and renderer design while retaining the production atomic path. Preserve historical measurement metadata and exclude generated Python bytecode.
191 lines
12 KiB
Plaintext
191 lines
12 KiB
Plaintext
/* Experimental pixel-owned reduction; intentionally not a production backend.
|
|
* Sharing this TU reuses the exact production cache indexing and row bounds. */
|
|
#include "../src/hip_psf.hip"
|
|
#include <omp.h>
|
|
#include <vector>
|
|
#include <algorithm>
|
|
#include <cstdlib>
|
|
#include <sys/resource.h>
|
|
|
|
struct TileTask { unsigned tile, first, count; };
|
|
__global__ static void tile_partial(const PsfCachedEvent *events, const unsigned *refs,
|
|
const TileTask *tasks, int tile_size, int tiles_x, int width, int height,
|
|
const float *weights, int phases, int radius, double max_radius, double *partial) {
|
|
const TileTask task=tasks[blockIdx.x];
|
|
const int local=threadIdx.x;
|
|
const int px=(task.tile%tiles_x)*tile_size+local%tile_size;
|
|
const int py=(task.tile/tiles_x)*tile_size+local/tile_size;
|
|
double r=0,g=0,b=0;
|
|
if(px<width && py<height) for(unsigned k=0;k<task.count;++k) {
|
|
const PsfCachedEvent e=events[refs[task.first+k]];
|
|
const int bx=(int)floor(e.x), by=(int)floor(e.y);
|
|
const int dx=px-bx, dy=py-by;
|
|
const int support=(int)fmin(ceil(e.support_radius),max_radius);
|
|
if(dy < -support || dy > support) continue;
|
|
const double fx=e.x-bx,fy=e.y-by;
|
|
int first,last;
|
|
if(!row_range(e.support_radius,fx,fy,support,dy,&first,&last) || dx<first || dx>last) continue;
|
|
const double ax=fx*phases,ay=fy*phases;
|
|
const int x0=(int)floor(ax),y0=(int)floor(ay);
|
|
const double tx=ax-x0,ty=ay-y0;
|
|
const double w00=weights[weight_index(phases,radius,x0,y0,dx,dy)];
|
|
const double w10=weights[weight_index(phases,radius,x0+1,y0,dx,dy)];
|
|
const double w01=weights[weight_index(phases,radius,x0,y0+1,dx,dy)];
|
|
const double w11=weights[weight_index(phases,radius,x0+1,y0+1,dx,dy)];
|
|
const double w=(1-ty)*((1-tx)*w00+tx*w10)+ty*((1-tx)*w01+tx*w11);
|
|
r+=e.color.r*e.flux*w;g+=e.color.g*e.flux*w;b+=e.color.b*e.flux*w;
|
|
}
|
|
const size_t p=3*((size_t)blockIdx.x*tile_size*tile_size+local);
|
|
partial[p]=r;partial[p+1]=g;partial[p+2]=b;
|
|
}
|
|
__global__ static void tile_merge(const unsigned *starts, int tile_size,int tiles_x,
|
|
int width,int height,const double *partial,double *hdr) {
|
|
const unsigned tile=blockIdx.x; const int local=threadIdx.x;
|
|
const int px=(tile%tiles_x)*tile_size+local%tile_size;
|
|
const int py=(tile/tiles_x)*tile_size+local/tile_size;
|
|
if(px>=width || py>=height) return;
|
|
double r=0,g=0,b=0;
|
|
for(unsigned k=starts[tile];k<starts[tile+1];++k) {
|
|
const size_t p=3*((size_t)k*tile_size*tile_size+local);
|
|
r+=partial[p];g+=partial[p+1];b+=partial[p+2];
|
|
}
|
|
const size_t p=3*((size_t)py*width+px);
|
|
hdr[p]+=r;hdr[p+1]+=g;hdr[p+2]+=b;
|
|
}
|
|
static void check(hipError_t e) {if(e!=hipSuccess){fprintf(stderr,"HIP: %s\n",hipGetErrorString(e));exit(1);}}
|
|
static double sync_time(HipPsfSink *s,double start) {check(hipStreamSynchronize(s->stream));return omp_get_wtime()-start;}
|
|
static size_t contributions(const std::vector<PsfCachedEvent>& events,int width,int height,const PsfKernelCache& cache) {
|
|
size_t n=0;
|
|
for(const auto&e:events) {
|
|
const int bx=(int)floor(e.x),by=(int)floor(e.y);
|
|
const double fx=e.x-bx,fy=e.y-by;
|
|
const int support=(int)fmin(ceil(e.support_radius),cache.max_radius_pixels);
|
|
for(int dy=-support;dy<=support;++dy) {
|
|
if(by+dy<0 || by+dy>=height) continue;
|
|
const double left=e.support_radius*e.support_radius-(dy+0.5-fy)*(dy+0.5-fy);
|
|
if(left<0) continue;
|
|
int lo=std::max(-support,(int)ceil(fx-0.5-sqrt(left)));
|
|
int hi=std::min(support,(int)floor(fx-0.5+sqrt(left)));
|
|
lo=std::max(lo,-bx);hi=std::min(hi,width-bx-1);
|
|
if(hi>=lo)n+=(size_t)(hi-lo+1);
|
|
}
|
|
}
|
|
return n;
|
|
}
|
|
int main(int argc,char **argv) {
|
|
if(argc<4 || argc>5) {puts("Usage: replay_psf INPUT.events EVENTS(1..65536) atomic|16|32 [disperse]\nExternal timeout <=45s; sequential processes only.");return 2;}
|
|
char *end=nullptr;const long requested=strtol(argv[2],&end,10);
|
|
if(!*argv[2] || *end || requested<1 || requested>65536) return 2;
|
|
const int tile=!strcmp(argv[3],"atomic") ? 0 : !strcmp(argv[3],"16") ? 16 : !strcmp(argv[3],"32") ? 32 : -1;
|
|
if(tile<0 || (argc==5 && strcmp(argv[4],"disperse")))return 2;
|
|
FILE *f=fopen(argv[1],"r");if(!f){perror(argv[1]);return 2;}
|
|
int width,height;size_t count;PointSpreadFunction psf;double tail;
|
|
if(fscanf(f,"PSFEVENTS1 %d %d %zu %lf %lf %lf",&width,&height,&count,&psf.fwhm_pixels,&psf.moffat_beta,&tail)!=6 ||
|
|
width<1 || height<1 || width>3840 || height>2160 || count>65536 || count<(size_t)requested ||
|
|
!std::isfinite(psf.fwhm_pixels) || psf.fwhm_pixels<=0 || psf.fwhm_pixels>6 ||
|
|
!std::isfinite(psf.moffat_beta) || psf.moffat_beta<2 || !std::isfinite(tail) || tail<1e-8 || tail>=1) return 2;
|
|
PsfKernelCache cache={};if(psf_kernel_cache_init(&cache,&psf,tail))return 1;
|
|
psf_kernel_cache_report_ready(&cache,stdout);
|
|
std::vector<PsfCachedEvent> events(requested);
|
|
for(auto &e:events) {
|
|
if(fscanf(f,"%lf %lf %lf %lf %lf %lf %lf",&e.x,&e.y,&e.color.r,&e.color.g,&e.color.b,&e.flux,&e.support_radius)!=7) return 2;
|
|
const double fields[]={e.x,e.y,e.color.r,e.color.g,e.color.b,e.flux,e.support_radius};
|
|
for(double v:fields)if(!std::isfinite(v))return 2;
|
|
if(fabs(e.x)>width+cache.max_radius_pixels || fabs(e.y)>height+cache.max_radius_pixels || e.support_radius<0 || e.support_radius>cache.max_radius_pixels)return 2;
|
|
}
|
|
fclose(f);
|
|
const size_t original_contributions=contributions(events,width,height,cache);
|
|
if(argc==5) {
|
|
unsigned state=12345;
|
|
for(auto&e:events) {
|
|
// Preserve fractional phase and clipping: move only fully interior events.
|
|
const int s=(int)ceil(e.support_radius)+1;
|
|
if(e.x<s || e.y<s || e.x>=width-s || e.y>=height-s)continue;
|
|
state=1664525u*state+1013904223u; e.x=s+state%(width-2*s)+(e.x-floor(e.x));
|
|
state=1664525u*state+1013904223u; e.y=s+state%(height-2*s)+(e.y-floor(e.y));
|
|
}
|
|
}
|
|
const size_t work=contributions(events,width,height,cache);
|
|
if(work!=original_contributions)return 1;
|
|
const size_t values=(size_t)width*height*3;
|
|
std::vector<double> cpu(values,0),gpu(values,0);
|
|
double start=omp_get_wtime();for(const auto&e:events)splat_prepared_cached_event(cpu.data(),width,height,&e,&cache);
|
|
const double cpu_time=omp_get_wtime()-start;
|
|
char message[256]={};HipPsfSink*s=nullptr;
|
|
start=omp_get_wtime();if(hip_psf_sink_create(&s,width,height,&cache,16384,message,sizeof message)){fprintf(stderr,"%s\n",message);return 1;}
|
|
const double create=omp_get_wtime()-start;
|
|
hipDeviceProp_t prop;check(hipGetDeviceProperties(&prop,0));
|
|
printf("device=%s arch=%s events=%zu frame=%dx%d mode=%s distribution=%s contributions=%zu CPU_reference=%.9f create=%.9f\n",prop.name,prop.gcnArchName,events.size(),width,height,argv[3],argc==5?"synthetic-phase-preserving-dispersal":"captured",work,cpu_time,create);fflush(stdout);
|
|
double bin=0,upload=0,kernel=0,merge=0,clear=0;size_t peak_scratch=0,ref_total=0,task_total=0;
|
|
unsigned *drefs=nullptr,*dstarts=nullptr;TileTask *dtasks=nullptr;double *partial=nullptr;
|
|
const double replay_start=omp_get_wtime();
|
|
if(tile) {
|
|
// Fixed hard bounds; no full-frame event list or per-chunk allocation on GPU.
|
|
check(hipMalloc(&drefs,32*1024*1024));check(hipMalloc(&dtasks,2*1024*1024));
|
|
check(hipMalloc(&dstarts,131072));check(hipMalloc(&partial,128*1024*1024));
|
|
}
|
|
for(size_t offset=0;offset<events.size();offset+=16384) {
|
|
const size_t n=std::min((size_t)16384,events.size()-offset);
|
|
if(!tile) {if(hip_psf_sink_submit(s,events.data()+offset,n,message,sizeof message)) {fprintf(stderr,"%s\n",message);return 1;}continue;}
|
|
start=omp_get_wtime();
|
|
const int nx=(width+tile-1)/tile,ny=(height+tile-1)/tile;
|
|
std::vector<std::vector<unsigned>> lists(nx*ny);
|
|
size_t refs_count=0;
|
|
for(size_t i=0;i<n;++i) {
|
|
const auto&e=events[offset+i];const int support=(int)fmin(ceil(e.support_radius),cache.max_radius_pixels);
|
|
const int bx=(int)floor(e.x),by=(int)floor(e.y);
|
|
const int x0=std::max(0,bx-support),x1=std::min(width-1,bx+support);
|
|
const int y0=std::max(0,by-support),y1=std::min(height-1,by+support);
|
|
if(x0>x1 || y0>y1)continue;
|
|
for(int y=y0/tile;y<=y1/tile;++y)for(int x=x0/tile;x<=x1/tile;++x) {
|
|
if(++refs_count>32*1024*1024/sizeof(unsigned)){fputs("reference limit exceeded\n",stderr);return 2;}
|
|
lists[y*nx+x].push_back((unsigned)i);
|
|
}
|
|
}
|
|
std::vector<unsigned> refs,starts(nx*ny+1);refs.reserve(refs_count);
|
|
std::vector<TileTask> tasks;
|
|
for(unsigned t=0;t<lists.size();++t) {
|
|
starts[t]=tasks.size();
|
|
for(size_t k=0;k<lists[t].size();k+=256) {
|
|
const unsigned length=std::min((size_t)256,lists[t].size()-k);
|
|
tasks.push_back({t,(unsigned)refs.size(),length});
|
|
refs.insert(refs.end(),lists[t].begin()+k,lists[t].begin()+k+length);
|
|
}
|
|
}
|
|
starts.back()=tasks.size();
|
|
const size_t bytes=tasks.size()*(size_t)tile*tile*3*sizeof(double);
|
|
if(bytes>128*1024*1024 || tasks.size()*sizeof(TileTask)>2*1024*1024 || starts.size()*sizeof(unsigned)>131072){fputs("partial/task limit exceeded\n",stderr);return 2;}
|
|
peak_scratch=std::max(peak_scratch,bytes+refs.size()*sizeof(unsigned)+tasks.size()*sizeof(TileTask)+starts.size()*sizeof(unsigned));
|
|
ref_total+=refs.size();task_total+=tasks.size();bin+=omp_get_wtime()-start;
|
|
start=omp_get_wtime();
|
|
check(hipMemcpyAsync(s->slots[0].device,events.data()+offset,n*sizeof(PsfCachedEvent),hipMemcpyHostToDevice,s->stream));
|
|
if(!refs.empty())check(hipMemcpyAsync(drefs,refs.data(),refs.size()*sizeof(unsigned),hipMemcpyHostToDevice,s->stream));
|
|
if(!tasks.empty())check(hipMemcpyAsync(dtasks,tasks.data(),tasks.size()*sizeof(TileTask),hipMemcpyHostToDevice,s->stream));
|
|
check(hipMemcpyAsync(dstarts,starts.data(),starts.size()*sizeof(unsigned),hipMemcpyHostToDevice,s->stream));
|
|
upload+=sync_time(s,start);
|
|
start=omp_get_wtime();
|
|
if(!tasks.empty())hipLaunchKernelGGL(tile_partial,dim3(tasks.size()),dim3(tile*tile),0,s->stream,s->slots[0].device,drefs,dtasks,tile,nx,width,height,s->weights,s->phase_resolution,s->radius_pixels,s->max_radius_pixels,partial);
|
|
check(hipGetLastError());kernel+=sync_time(s,start);
|
|
start=omp_get_wtime();hipLaunchKernelGGL(tile_merge,dim3(nx*ny),dim3(tile*tile),0,s->stream,dstarts,tile,nx,width,height,partial,s->hdr);
|
|
check(hipGetLastError());merge+=sync_time(s,start);
|
|
// Every partial pixel is written by exactly one thread: no memset required.
|
|
}
|
|
if(hip_psf_sink_finish(s,gpu.data(),message,sizeof message)){fprintf(stderr,"%s\n",message);return 1;}
|
|
if(tile) {start=omp_get_wtime();check(hipFree(drefs));check(hipFree(dtasks));check(hipFree(dstarts));check(hipFree(partial));clear=omp_get_wtime()-start;}
|
|
const double wall=omp_get_wtime()-replay_start;
|
|
HipPsfTiming timing={};hip_psf_sink_get_timing(s,&timing);
|
|
if(!tile){upload=timing.upload_seconds;kernel=timing.kernel_seconds;}
|
|
double max_abs=0,max_rel=0;long double sums[3]={},diffs[3]={};bool pass=true;
|
|
for(size_t i=0;i<values;++i) {
|
|
if(!std::isfinite(cpu[i]) || !std::isfinite(gpu[i]))pass=false;
|
|
const double d=fabs(gpu[i]-cpu[i]);max_abs=fmax(max_abs,d);
|
|
if(cpu[i]!=0)max_rel=fmax(max_rel,d/fabs(cpu[i]));
|
|
sums[i%3]+=cpu[i];diffs[i%3]+=(long double)gpu[i]-cpu[i];
|
|
}
|
|
long double ysum=.2126L*sums[0]+.7152L*sums[1]+.0722L*sums[2],ydiff=.2126L*diffs[0]+.7152L*diffs[1]+.0722L*diffs[2];
|
|
struct rusage usage;getrusage(RUSAGE_SELF,&usage);
|
|
printf("replay_wall=%.9f bin=%.9f upload=%.9f accumulation=%.9f merge=%.9f download=%.9f cleanup=%.9f events_s=%.3f contributions_s=%.3f refs=%zu tasks=%zu scratch_used_peak=%zu scratch_device_reserved=%u RSS_KiB=%ld max_abs=%.17g max_rel=%.17g flux_R=%.17Lg flux_G=%.17Lg flux_B=%.17Lg flux_Y=%.17Lg\n",wall,bin,upload,kernel,merge,timing.download_seconds,clear,events.size()/wall,work/wall,ref_total,task_total,peak_scratch,tile?162*1024*1024+131072:0,usage.ru_maxrss,max_abs,max_rel,sums[0]?diffs[0]/sums[0]:0,sums[1]?diffs[1]/sums[1]:0,sums[2]?diffs[2]/sums[2]:0,ysum?ydiff/ysum:0);
|
|
hip_psf_sink_destroy(s);psf_kernel_cache_destroy(&cache);
|
|
return pass && max_abs<1e-10 && max_rel<1e-10 ? 0 : 1;
|
|
}
|