Skip to content

Commit 6a43448

Browse files
falseywinchnetfalseywinchnet
authored
Merge pull request #67 from falseywinchnet/jmprpn-main/update-mag+phase-output-process
Delay mag-phase conversion and wire atan benchmark
2 parents 63e927a + 225a8a9 commit 6a43448

3 files changed

Lines changed: 84 additions & 29 deletions

File tree

‎Makefile‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -57,6 +57,7 @@ PC_FILE := $(BUILD_DIR)/$(LIB_NAME).pc
5757
BENCH := $(BUILD_DIR)/examples/benchmark
5858
BODFT_BENCH := $(BUILD_DIR)/examples/bodft_benchmark
5959
APPLE_BENCH := $(BUILD_DIR)/examples/apple_benchmark
60+
ATAN2_BENCH := $(BUILD_DIR)/examples/atan2_benchmark
6061
LOCALITY_PROBE := $(BUILD_DIR)/examples/locality_probe
6162
C_DEMO := $(BUILD_DIR)/examples/c_api_demo
6263
CPP_DEMO := $(BUILD_DIR)/examples/cpp_api_demo
@@ -109,7 +110,7 @@ $(PC_FILE): pkgconfig/bfft.pc.in | $(BUILD_DIR)
109110

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

112-
benchmarks: $(BENCH) $(BODFT_BENCH)
113+
benchmarks: $(BENCH) $(BODFT_BENCH) $(ATAN2_BENCH)
113114

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

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

134136
$(APPLE_BENCH): examples/apple_benchmark.cpp include/bfft/bfft.hpp $(STATIC_LIB) | $(BUILD_DIR)
135137
$(CXX) $(CPPFLAGS) $(INCLUDES) $(CXXFLAGS) $< $(STATIC_LIB) $(LDLIBS) $(DL_LIBS) $(ACCELERATE_LIBS) -o $@

‎examples/atan2_benchmark.cpp‎

Lines changed: 22 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -62,6 +62,26 @@ double time_loop(const std::vector<sample>& samples, Func func, double& sink) {
6262
return std::chrono::duration<double>(stop - start).count();
6363
}
6464

65+
double time_vec5_pair_loop(const std::vector<sample>& samples, double& sink) {
66+
const auto start = std::chrono::steady_clock::now();
67+
double total = 0.0;
68+
std::size_t i = 0;
69+
for (; i + 1 < samples.size(); i += 2) {
70+
double phase0;
71+
double phase1;
72+
bruun::bruun_phase_atan2_mag_pair(samples[i].y, samples[i].x, samples[i].mag,
73+
samples[i + 1].y, samples[i + 1].x, samples[i + 1].mag,
74+
&phase0, &phase1);
75+
total += phase0 + phase1;
76+
}
77+
for (; i < samples.size(); ++i) {
78+
total += bruun::bruun_phase_atan2_mag(samples[i].y, samples[i].x, samples[i].mag);
79+
}
80+
const auto stop = std::chrono::steady_clock::now();
81+
sink += total;
82+
return std::chrono::duration<double>(stop - start).count();
83+
}
84+
6585
} // namespace
6686

6787
int main(int argc, char** argv) {
@@ -82,9 +102,7 @@ int main(int argc, char** argv) {
82102
bruun::bruun_phase_first_octant_tree32_degree7);
83103
}, sink);
84104

85-
const double bfft_vec_double_seconds = time_loop(samples, [](const sample& value) {
86-
return bruun::bruun_phase_atan2_mag(value.y, value.x, value.mag);
87-
}, sink);
105+
const double bfft_vec_double_seconds = time_vec5_pair_loop(samples, sink);
88106

89107
const double std_float_seconds = time_loop(samples, [](const sample& value) {
90108
const float y = static_cast<float>(value.y);
@@ -121,7 +139,7 @@ int main(int argc, char** argv) {
121139
std::printf("samples=%zu sink=%.17g\n", count, sink);
122140
std::printf("std::atan2 double: %.6f s, %.3f Mphase/s\n", std_double_seconds, scale / std_double_seconds);
123141
std::printf("bfft tree32 double: %.6f s, %.3f Mphase/s\n", bfft_tree_double_seconds, scale / bfft_tree_double_seconds);
124-
std::printf("bfft vec5 double : %.6f s, %.3f Mphase/s\n", bfft_vec_double_seconds, scale / bfft_vec_double_seconds);
142+
std::printf("bfft vec5 pair : %.6f s, %.3f Mphase/s\n", bfft_vec_double_seconds, scale / bfft_vec_double_seconds);
125143
std::printf("std::atan2 float : %.6f s, %.3f Mphase/s\n", std_float_seconds, scale / std_float_seconds);
126144
std::printf("bfft tree32 float : %.6f s, %.3f Mphase/s\n", bfft_tree_float_seconds, scale / bfft_tree_float_seconds);
127145
std::printf("bfft vec5 float : %.6f s, %.3f Mphase/s\n", bfft_vec_float_seconds, scale / bfft_vec_float_seconds);

‎src/detail/bruun_kernel.hpp‎

Lines changed: 58 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -2736,15 +2736,23 @@ class RFFT {
27362736
}
27372737

27382738
const int* RESTRICT kin = KINV.data();
2739+
for (int k = 1; k < N / 2; ++k) {
2740+
const int m = kin[k];
2741+
output[k].re = work[2*m];
2742+
output[k].im = -work[2*m + 1];
2743+
}
2744+
2745+
convert_standard_complex_to_mag_phase(output);
2746+
}
2747+
2748+
void convert_standard_complex_to_mag_phase(complex_t* RESTRICT output) const {
27392749
int k = 1;
27402750
#if BRUUN_LEVEL >= 1
27412751
for (; k + 1 < N / 2; k += 2) {
2742-
const int m0 = kin[k];
2743-
const int m1 = kin[k + 1];
2744-
const double re0 = work[2*m0];
2745-
const double im0 = -work[2*m0 + 1];
2746-
const double re1 = work[2*m1];
2747-
const double im1 = -work[2*m1 + 1];
2752+
const double re0 = output[k].re;
2753+
const double im0 = output[k].im;
2754+
const double re1 = output[k + 1].re;
2755+
const double im1 = output[k + 1].im;
27482756
const double mag0 = std::sqrt(re0 * re0 + im0 * im0);
27492757
const double mag1 = std::sqrt(re1 * re1 + im1 * im1);
27502758
double phase0;
@@ -2757,14 +2765,41 @@ class RFFT {
27572765
}
27582766
#endif
27592767
for (; k < N / 2; ++k) {
2760-
const int m = kin[k];
2761-
const double re = work[2*m];
2762-
const double im = -work[2*m + 1];
2768+
const double re = output[k].re;
2769+
const double im = output[k].im;
27632770
const double mag = std::sqrt(re * re + im * im);
2764-
double phase = bruun_phase_atan2_mag(im, re, mag);
2765-
if (phase < 0.0) {
2766-
phase += 2.0 * M_PI;
2767-
}
2771+
const double phase = bruun_phase_atan2_mag(im, re, mag);
2772+
output[k].re = mag;
2773+
output[k].im = phase;
2774+
}
2775+
}
2776+
2777+
void convert_standard_complex_to_mag_phase_f32(complex_f32_t* RESTRICT output) const {
2778+
int k = 1;
2779+
#if BRUUN_LEVEL >= 1
2780+
for (; k + 1 < N / 2; k += 2) {
2781+
const float re0 = output[k].re;
2782+
const float im0 = output[k].im;
2783+
const float re1 = output[k + 1].re;
2784+
const float im1 = output[k + 1].im;
2785+
const float mag0 = std::sqrt(re0 * re0 + im0 * im0);
2786+
const float mag1 = std::sqrt(re1 * re1 + im1 * im1);
2787+
double phase0;
2788+
double phase1;
2789+
bruun_phase_atan2_mag_pair(static_cast<double>(im0), static_cast<double>(re0), static_cast<double>(mag0),
2790+
static_cast<double>(im1), static_cast<double>(re1), static_cast<double>(mag1),
2791+
&phase0, &phase1);
2792+
output[k].re = mag0;
2793+
output[k].im = static_cast<float>(phase0);
2794+
output[k + 1].re = mag1;
2795+
output[k + 1].im = static_cast<float>(phase1);
2796+
}
2797+
#endif
2798+
for (; k < N / 2; ++k) {
2799+
const float re = output[k].re;
2800+
const float im = output[k].im;
2801+
const float mag = std::sqrt(re * re + im * im);
2802+
const float phase = bruun_phase_atan2_mag_f32(im, re, mag);
27682803
output[k].re = mag;
27692804
output[k].im = phase;
27702805
}
@@ -2812,18 +2847,18 @@ class RFFT {
28122847
}
28132848

28142849
const int* RESTRICT kin = KINV.data();
2815-
for (int k = 1; k < N / 2; ++k) {
2850+
int k = 1;
2851+
for (; k + 3 < N / 2; k += 4) {
2852+
pack4_residues_to_complex_f32(output + k, work,
2853+
kin[k], kin[k + 1], kin[k + 2], kin[k + 3]);
2854+
}
2855+
for (; k < N / 2; ++k) {
28162856
const int m = kin[k];
2817-
const float re = work[2*m];
2818-
const float im = -work[2*m + 1];
2819-
const float mag = std::sqrt(re * re + im * im);
2820-
float phase = bruun_phase_atan2_mag_f32(im, re, mag);
2821-
if (phase < 0.0f) {
2822-
phase += 2.0f * static_cast<float>(M_PI);
2823-
}
2824-
output[k].re = mag;
2825-
output[k].im = phase;
2857+
output[k].re = work[2*m];
2858+
output[k].im = -work[2*m + 1];
28262859
}
2860+
2861+
convert_standard_complex_to_mag_phase_f32(output);
28272862
}
28282863

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

0 commit comments

Comments
 (0)