diff --git a/psf_lookup_optimization_plan.md b/psf_lookup_optimization_plan.md index 028410e..ff6e43d 100644 --- a/psf_lookup_optimization_plan.md +++ b/psf_lookup_optimization_plan.md @@ -116,7 +116,8 @@ case. The direct evaluator remains a selectable regression reference. Implemented as one immutable `PsfKernelCache` per invocation, built after CLI parsing and passed read-only through the existing frame splat workers. The cache uses 64 by 64 phase intervals (65 nodes on each axis), 4-point -Gauss-Legendre pixel-area quadrature, and bilinear phase interpolation. The +Gauss-Legendre pixel-area quadrature, bilinear phase interpolation, and a +static OpenMP phase-plane construction loop. The `--psf-direct` regression mode uses an independent 8-point quadrature reference; events whose requested support exceeds the cache use that path automatically. diff --git a/src/optics.c b/src/optics.c index 6bcb087..0465d32 100644 --- a/src/optics.c +++ b/src/optics.c @@ -2,11 +2,11 @@ #include #include +#include #include #include #include #include -#include #ifdef ENABLE_PNG #include @@ -213,7 +213,7 @@ int psf_kernel_cache_init(PsfKernelCache *cache, const PointSpreadFunction *psf) float *weights = malloc(count * sizeof *weights); if (weights == NULL) return -1; - const clock_t start = clock(); + const double start = omp_get_wtime(); PsfKernelCache building = {.fwhm_pixels = psf->fwhm_pixels, .moffat_beta = psf->moffat_beta, .alpha_pixels = alpha, @@ -221,6 +221,11 @@ int psf_kernel_cache_init(PsfKernelCache *cache, const PointSpreadFunction *psf) .weights = weights, .phase_resolution = PSF_PHASE_RESOLUTION, .radius_pixels = radius}; + int invalid_kernel = 0; + /* Each sub-pixel phase owns a disjoint kernel plane. Parallelizing this + * regular construction needs no locks, and the finished kernel remains + * immutable to splat workers. */ +#pragma omp parallel for collapse(2) schedule(static) for (int phase_y = 0; phase_y <= PSF_PHASE_RESOLUTION; ++phase_y) for (int phase_x = 0; phase_x <= PSF_PHASE_RESOLUTION; ++phase_x) { const double star_x = (double)phase_x / PSF_PHASE_RESOLUTION; @@ -234,10 +239,12 @@ int psf_kernel_cache_init(PsfKernelCache *cache, const PointSpreadFunction *psf) alpha, psf->moffat_beta, offset_x, offset_y, star_x, star_y); weights[index] = (float)weight; sum += weight; - } + } if (!(sum > 0.0) || !isfinite(sum)) { - free(weights); - return -1; + /* An error cannot return from an OpenMP structured block. */ +#pragma omp atomic write + invalid_kernel = 1; + continue; } for (int offset_y = -radius; offset_y <= radius; ++offset_y) for (int offset_x = -radius; offset_x <= radius; ++offset_x) { @@ -246,7 +253,11 @@ int psf_kernel_cache_init(PsfKernelCache *cache, const PointSpreadFunction *psf) weights[index] = (float)(weights[index] / sum); } } - building.build_seconds = (double)(clock() - start) / CLOCKS_PER_SEC; + if (invalid_kernel) { + free(weights); + return -1; + } + building.build_seconds = omp_get_wtime() - start; building.ready = 1; *cache = building; return 0;