Commit f5373a50 authored by Leo Pound Singer's avatar Leo Pound Singer
Browse files

Reduce the number of log() and exp() calls in inner loop

For large numbers of terms, logaddexp() performs about twice as many
transcendental function evaluations as is necessary.
This patch results in a speedup of about 25%.
Original: 41634d246522aefd772d5f6b19d13277476b9c07
parent 8947e0d1
......@@ -832,7 +832,11 @@ bayestar_pixel *bayestar_sky_map_toa_phoa_snr(
for (unsigned int iu = 0; iu < nglfixed; iu++)
{
double u, weight;
double accum2[3] = {-INFINITY, -INFINITY, -INFINITY};
double accum2[nsamples][3];
for (long isample = 0; isample < 3; isample ++)
for (unsigned char k = 0; k < 3; k ++)
accum2[isample][k] = -INFINITY;
{
/* Look up Gauss-Legendre abscissa and weight. */
int ret = gsl_integration_glfixed_point(
......@@ -874,16 +878,22 @@ bayestar_pixel *bayestar_sky_map_toa_phoa_snr(
for (unsigned char k = 0; k < 3; k ++)
{
double result = log_radial_integrator_eval(
accum2[isample][k] = log_radial_integrator_eval(
integrators[k], p, b, log_p, log_b);
accum2[k] = logaddexp(accum2[k], result);
}
}
double max_accum2[3] = {-INFINITY, -INFINITY, -INFINITY};
for (long isample = 0; isample < 3; isample ++)
for (unsigned char k = 0; k < 3; k ++)
if (accum2[isample][k] > max_accum2[k])
max_accum2[k] = accum2[isample][k];
double sum_accum2[3] = {0, 0, 0};
for (long isample = 0; isample < 3; isample ++)
for (unsigned char k = 0; k < 3; k ++)
sum_accum2[k] += exp(accum2[isample][k] - max_accum2[k]);
for (unsigned char k = 0; k < 3; k ++)
{
accum1[k] = logaddexp(accum1[k], accum2[k] + log(weight));
}
accum1[k] = logaddexp(accum1[k], log(sum_accum2[k] * weight));
}
for (unsigned char k = 0; k < 3; k ++)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment