Skip to content

spec_avg.c: make sure SFT real/imag components are REAL8 before multiplying

Karl Wette requested to merge ANU-CGA/lalsuite:icc-test into master

Description

Fixes the failing testFscan.sh under ICC.

The problem was this line in lalapps/src/pulsar/Fscan/spec_avg.c:

avg += 2.0*(crealf(sft_vect->data[j].data->data[i-k])*crealf(sft_vect->data[j].data->data[i-k])
          + cimagf(sft_vect->data[j].data->data[i-k])*cimagf(sft_vect->data[j].data->data[i-k]))/(REAL8)timebaseline;

The SFT data sft_vect->data[j].data->data[i-k] is a COMPLEX8, so in terms of types this line is

REAL8 += 2.0*(REAL4*REAL4 + REAL4*REAL4)/REAL8;

The problem is that ICC isn't upconverting REAL4 values to REAL8 when assigning to REAL8 (whereas GCC evidently does), which leads to loss of precision. In particular, SFT data is of the order of h(t) ~ 1e-21, and (1e-21)^2 is too small to hold in REAL4 precision, so squaring the SFT data in REAL4 results in zeros being added to avg.

The solution is to assign the real and imaginary parts of the SFT data explicitly to REAL8 first, then do the multiplication:

const REAL8 re = (REAL8)crealf(sft_vect->data[j].data->data[i-k]);
const REAL8 im = (REAL8)cimagf(sft_vect->data[j].data->data[i-k]);
avg += 2.0*(re*re + im*im)/(REAL8)timebaseline;

(Which has the bonus of being a far more compact expression.)

To be safe I made the same change for the calculation involving timeavg->data[i] later in spec_avg.c

Closes #460 (closed)

API Changes and Justification

Backwards Compatible Changes

  • This change introduces no API changes
  • This change adds new API calls

Backwards Incompatible Changes

  • This change modifies an existing API
  • This change removes an existing API

Review Status

cc @evan-goetz @adam-mercer

Merge request reports