Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 4 additions & 2 deletions Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,7 @@ PC_FILE := $(BUILD_DIR)/$(LIB_NAME).pc
BENCH := $(BUILD_DIR)/examples/benchmark
BODFT_BENCH := $(BUILD_DIR)/examples/bodft_benchmark
APPLE_BENCH := $(BUILD_DIR)/examples/apple_benchmark
ATAN2_BENCH := $(BUILD_DIR)/examples/atan2_benchmark
LOCALITY_PROBE := $(BUILD_DIR)/examples/locality_probe
C_DEMO := $(BUILD_DIR)/examples/c_api_demo
CPP_DEMO := $(BUILD_DIR)/examples/cpp_api_demo
Expand Down Expand Up @@ -109,7 +110,7 @@ $(PC_FILE): pkgconfig/bfft.pc.in | $(BUILD_DIR)

examples: $(BENCH) $(BODFT_BENCH) $(APPLE_EXAMPLES) $(LOCALITY_PROBE) $(C_DEMO) $(CPP_DEMO)

benchmarks: $(BENCH) $(BODFT_BENCH)
benchmarks: $(BENCH) $(BODFT_BENCH) $(ATAN2_BENCH)

asm-check: $(ASM_OUTPUTS)
@if [ -z "$(ASM_OUTPUTS)" ]; then echo "No x86 assembly variants supported by $(CXX)."; fi
Expand All @@ -129,7 +130,8 @@ $(BENCH): examples/benchmark.cpp include/bfft/bfft.hpp $(STATIC_LIB) | $(BUILD_D
$(BODFT_BENCH): examples/bodft_benchmark.cpp src/detail/bodft_kernel.hpp src/detail/bruun_kernel.hpp | $(BUILD_DIR)
$(CXX) $(CPPFLAGS) $(INCLUDES) $(CXXFLAGS) $(AUTO_SIMD_FLAGS) $< $(LDLIBS) -o $@


$(ATAN2_BENCH): examples/atan2_benchmark.cpp src/detail/bruun_kernel.hpp | $(BUILD_DIR)
$(CXX) $(CPPFLAGS) $(INCLUDES) $(CXXFLAGS) $(AUTO_SIMD_FLAGS) $< $(LDLIBS) -o $@

$(APPLE_BENCH): examples/apple_benchmark.cpp include/bfft/bfft.hpp $(STATIC_LIB) | $(BUILD_DIR)
$(CXX) $(CPPFLAGS) $(INCLUDES) $(CXXFLAGS) $< $(STATIC_LIB) $(LDLIBS) $(DL_LIBS) $(ACCELERATE_LIBS) -o $@
Expand Down
26 changes: 22 additions & 4 deletions examples/atan2_benchmark.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,26 @@ double time_loop(const std::vector<sample>& samples, Func func, double& sink) {
return std::chrono::duration<double>(stop - start).count();
}

double time_vec5_pair_loop(const std::vector<sample>& samples, double& sink) {
const auto start = std::chrono::steady_clock::now();
double total = 0.0;
std::size_t i = 0;
for (; i + 1 < samples.size(); i += 2) {
double phase0;
double phase1;
bruun::bruun_phase_atan2_mag_pair(samples[i].y, samples[i].x, samples[i].mag,
samples[i + 1].y, samples[i + 1].x, samples[i + 1].mag,
&phase0, &phase1);
total += phase0 + phase1;
}
for (; i < samples.size(); ++i) {
total += bruun::bruun_phase_atan2_mag(samples[i].y, samples[i].x, samples[i].mag);
}
const auto stop = std::chrono::steady_clock::now();
sink += total;
return std::chrono::duration<double>(stop - start).count();
}

} // namespace

int main(int argc, char** argv) {
Expand All @@ -82,9 +102,7 @@ int main(int argc, char** argv) {
bruun::bruun_phase_first_octant_tree32_degree7);
}, sink);

const double bfft_vec_double_seconds = time_loop(samples, [](const sample& value) {
return bruun::bruun_phase_atan2_mag(value.y, value.x, value.mag);
}, sink);
const double bfft_vec_double_seconds = time_vec5_pair_loop(samples, sink);

const double std_float_seconds = time_loop(samples, [](const sample& value) {
const float y = static_cast<float>(value.y);
Expand Down Expand Up @@ -121,7 +139,7 @@ int main(int argc, char** argv) {
std::printf("samples=%zu sink=%.17g\n", count, sink);
std::printf("std::atan2 double: %.6f s, %.3f Mphase/s\n", std_double_seconds, scale / std_double_seconds);
std::printf("bfft tree32 double: %.6f s, %.3f Mphase/s\n", bfft_tree_double_seconds, scale / bfft_tree_double_seconds);
std::printf("bfft vec5 double : %.6f s, %.3f Mphase/s\n", bfft_vec_double_seconds, scale / bfft_vec_double_seconds);
std::printf("bfft vec5 pair : %.6f s, %.3f Mphase/s\n", bfft_vec_double_seconds, scale / bfft_vec_double_seconds);
std::printf("std::atan2 float : %.6f s, %.3f Mphase/s\n", std_float_seconds, scale / std_float_seconds);
std::printf("bfft tree32 float : %.6f s, %.3f Mphase/s\n", bfft_tree_float_seconds, scale / bfft_tree_float_seconds);
std::printf("bfft vec5 float : %.6f s, %.3f Mphase/s\n", bfft_vec_float_seconds, scale / bfft_vec_float_seconds);
Expand Down
81 changes: 58 additions & 23 deletions src/detail/bruun_kernel.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -2736,15 +2736,23 @@ class RFFT {
}

const int* RESTRICT kin = KINV.data();
for (int k = 1; k < N / 2; ++k) {
const int m = kin[k];
output[k].re = work[2*m];
output[k].im = -work[2*m + 1];
}

convert_standard_complex_to_mag_phase(output);
}

void convert_standard_complex_to_mag_phase(complex_t* RESTRICT output) const {
int k = 1;
#if BRUUN_LEVEL >= 1
for (; k + 1 < N / 2; k += 2) {
const int m0 = kin[k];
const int m1 = kin[k + 1];
const double re0 = work[2*m0];
const double im0 = -work[2*m0 + 1];
const double re1 = work[2*m1];
const double im1 = -work[2*m1 + 1];
const double re0 = output[k].re;
const double im0 = output[k].im;
const double re1 = output[k + 1].re;
const double im1 = output[k + 1].im;
const double mag0 = std::sqrt(re0 * re0 + im0 * im0);
const double mag1 = std::sqrt(re1 * re1 + im1 * im1);
double phase0;
Expand All @@ -2757,14 +2765,41 @@ class RFFT {
}
#endif
for (; k < N / 2; ++k) {
const int m = kin[k];
const double re = work[2*m];
const double im = -work[2*m + 1];
const double re = output[k].re;
const double im = output[k].im;
const double mag = std::sqrt(re * re + im * im);
double phase = bruun_phase_atan2_mag(im, re, mag);
if (phase < 0.0) {
phase += 2.0 * M_PI;
}
const double phase = bruun_phase_atan2_mag(im, re, mag);
output[k].re = mag;
output[k].im = phase;
}
}

void convert_standard_complex_to_mag_phase_f32(complex_f32_t* RESTRICT output) const {
int k = 1;
#if BRUUN_LEVEL >= 1
for (; k + 1 < N / 2; k += 2) {
const float re0 = output[k].re;
const float im0 = output[k].im;
const float re1 = output[k + 1].re;
const float im1 = output[k + 1].im;
const float mag0 = std::sqrt(re0 * re0 + im0 * im0);
const float mag1 = std::sqrt(re1 * re1 + im1 * im1);
double phase0;
double phase1;
bruun_phase_atan2_mag_pair(static_cast<double>(im0), static_cast<double>(re0), static_cast<double>(mag0),
static_cast<double>(im1), static_cast<double>(re1), static_cast<double>(mag1),
&phase0, &phase1);
output[k].re = mag0;
output[k].im = static_cast<float>(phase0);
output[k + 1].re = mag1;
output[k + 1].im = static_cast<float>(phase1);
}
#endif
for (; k < N / 2; ++k) {
const float re = output[k].re;
const float im = output[k].im;
const float mag = std::sqrt(re * re + im * im);
const float phase = bruun_phase_atan2_mag_f32(im, re, mag);
output[k].re = mag;
output[k].im = phase;
}
Expand Down Expand Up @@ -2812,18 +2847,18 @@ class RFFT {
}

const int* RESTRICT kin = KINV.data();
for (int k = 1; k < N / 2; ++k) {
int k = 1;
for (; k + 3 < N / 2; k += 4) {
pack4_residues_to_complex_f32(output + k, work,
kin[k], kin[k + 1], kin[k + 2], kin[k + 3]);
}
for (; k < N / 2; ++k) {
const int m = kin[k];
const float re = work[2*m];
const float im = -work[2*m + 1];
const float mag = std::sqrt(re * re + im * im);
float phase = bruun_phase_atan2_mag_f32(im, re, mag);
if (phase < 0.0f) {
phase += 2.0f * static_cast<float>(M_PI);
}
output[k].re = mag;
output[k].im = phase;
output[k].re = work[2*m];
output[k].im = -work[2*m + 1];
}

convert_standard_complex_to_mag_phase_f32(output);
}

// Standard FFTW-like real -> complex interface, using caller-provided native scratch.
Expand Down
Loading