From 7887b1044b0de7d6c568336f52fd0eab5e7f643d Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Fri, 12 Jun 2020 15:56:45 -0400 Subject: [PATCH] ensure highly factorizable MPB sizes (#1244) --- src/mpb.cpp | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/mpb.cpp b/src/mpb.cpp index cc071899e..e4fa3875d 100644 --- a/src/mpb.cpp +++ b/src/mpb.cpp @@ -233,6 +233,20 @@ static double dot_product(const mpb_real a[3], const mpb_real b[3]) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; } +// return the next number factorizable into powers of 2,3,5,7 +// for efficient FFTs, similar to optimize-grid-size! in MPB: +static int nextpow2357(int n) { + while (1) { + int m = n; + while (m % 3 == 0) m /= 3; + while (m % 5 == 0) m /= 5; + while (m % 7 == 0) m /= 7; + if ((m & (m - 1)) == 0) // if m is a power of 2 + return n; + n += 1; + } +} + /****************************************************************/ /* call MPB to get the band_numth eigenmode at freq frequency. */ /* */ @@ -332,6 +346,7 @@ void *fields::get_eigenmode(double frequency, direction d, const volume where, c for (int i = 0; i < 3; ++i) { n[i] = int(resolution * s[i] + 0.5); if (n[i] == 0) n[i] = 1; + n[i] = nextpow2357(n[i]); if (s[i] != 0) R[i][i] = s[i]; else {