Mercurial > hg > sv-dependency-builds

annotate src/fftw-3.3.5/dft/bluestein.c @ 84:08ae793730bd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add null config files
author Chris Cannam
date Mon, 02 Mar 2020 14:03:47 +0000
parents 2cd0e3b3e1fd
children
rev   line source
Chris@42 1 /*
Chris@42 2 * Copyright (c) 2003, 2007-14 Matteo Frigo
Chris@42 3 * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
Chris@42 4 *
Chris@42 5 * This program is free software; you can redistribute it and/or modify
Chris@42 6 * it under the terms of the GNU General Public License as published by
Chris@42 7 * the Free Software Foundation; either version 2 of the License, or
Chris@42 8 * (at your option) any later version.
Chris@42 9 *
Chris@42 10 * This program is distributed in the hope that it will be useful,
Chris@42 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@42 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
Chris@42 13 * GNU General Public License for more details.
Chris@42 14 *
Chris@42 15 * You should have received a copy of the GNU General Public License
Chris@42 16 * along with this program; if not, write to the Free Software
Chris@42 17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
Chris@42 18 *
Chris@42 19 */
Chris@42 20
Chris@42 21 #include "dft.h"
Chris@42 22
Chris@42 23 typedef struct {
Chris@42 24 solver super;
Chris@42 25 } S;
Chris@42 26
Chris@42 27 typedef struct {
Chris@42 28 plan_dft super;
Chris@42 29 INT n; /* problem size */
Chris@42 30 INT nb; /* size of convolution */
Chris@42 31 R *w; /* lambda k . exp(2*pi*i*k^2/(2*n)) */
Chris@42 32 R *W; /* DFT(w) */
Chris@42 33 plan *cldf;
Chris@42 34 INT is, os;
Chris@42 35 } P;
Chris@42 36
Chris@42 37 static void bluestein_sequence(enum wakefulness wakefulness, INT n, R *w)
Chris@42 38 {
Chris@42 39 INT k, ksq, n2 = 2 * n;
Chris@42 40 triggen *t = X(mktriggen)(wakefulness, n2);
Chris@42 41
Chris@42 42 ksq = 0;
Chris@42 43 for (k = 0; k < n; ++k) {
Chris@42 44 t->cexp(t, ksq, w+2*k);
Chris@42 45 /* careful with overflow */
Chris@42 46 ksq += 2*k + 1; while (ksq > n2) ksq -= n2;
Chris@42 47 }
Chris@42 48
Chris@42 49 X(triggen_destroy)(t);
Chris@42 50 }
Chris@42 51
Chris@42 52 static void mktwiddle(enum wakefulness wakefulness, P *p)
Chris@42 53 {
Chris@42 54 INT i;
Chris@42 55 INT n = p->n, nb = p->nb;
Chris@42 56 R *w, *W;
Chris@42 57 E nbf = (E)nb;
Chris@42 58
Chris@42 59 p->w = w = (R *) MALLOC(2 * n * sizeof(R), TWIDDLES);
Chris@42 60 p->W = W = (R *) MALLOC(2 * nb * sizeof(R), TWIDDLES);
Chris@42 61
Chris@42 62 bluestein_sequence(wakefulness, n, w);
Chris@42 63
Chris@42 64 for (i = 0; i < nb; ++i)
Chris@42 65 W[2*i] = W[2*i+1] = K(0.0);
Chris@42 66
Chris@42 67 W[0] = w[0] / nbf;
Chris@42 68 W[1] = w[1] / nbf;
Chris@42 69
Chris@42 70 for (i = 1; i < n; ++i) {
Chris@42 71 W[2*i] = W[2*(nb-i)] = w[2*i] / nbf;
Chris@42 72 W[2*i+1] = W[2*(nb-i)+1] = w[2*i+1] / nbf;
Chris@42 73 }
Chris@42 74
Chris@42 75 {
Chris@42 76 plan_dft *cldf = (plan_dft *)p->cldf;
Chris@42 77 /* cldf must be awake */
Chris@42 78 cldf->apply(p->cldf, W, W+1, W, W+1);
Chris@42 79 }
Chris@42 80 }
Chris@42 81
Chris@42 82 static void apply(const plan *ego_, R *ri, R *ii, R *ro, R *io)
Chris@42 83 {
Chris@42 84 const P *ego = (const P *) ego_;
Chris@42 85 INT i, n = ego->n, nb = ego->nb, is = ego->is, os = ego->os;
Chris@42 86 R *w = ego->w, *W = ego->W;
Chris@42 87 R *b = (R *) MALLOC(2 * nb * sizeof(R), BUFFERS);
Chris@42 88
Chris@42 89 /* multiply input by conjugate bluestein sequence */
Chris@42 90 for (i = 0; i < n; ++i) {
Chris@42 91 E xr = ri[i*is], xi = ii[i*is];
Chris@42 92 E wr = w[2*i], wi = w[2*i+1];
Chris@42 93 b[2*i] = xr * wr + xi * wi;
Chris@42 94 b[2*i+1] = xi * wr - xr * wi;
Chris@42 95 }
Chris@42 96
Chris@42 97 for (; i < nb; ++i) b[2*i] = b[2*i+1] = K(0.0);
Chris@42 98
Chris@42 99 /* convolution: FFT */
Chris@42 100 {
Chris@42 101 plan_dft *cldf = (plan_dft *)ego->cldf;
Chris@42 102 cldf->apply(ego->cldf, b, b+1, b, b+1);
Chris@42 103 }
Chris@42 104
Chris@42 105 /* convolution: pointwise multiplication */
Chris@42 106 for (i = 0; i < nb; ++i) {
Chris@42 107 E xr = b[2*i], xi = b[2*i+1];
Chris@42 108 E wr = W[2*i], wi = W[2*i+1];
Chris@42 109 b[2*i] = xi * wr + xr * wi;
Chris@42 110 b[2*i+1] = xr * wr - xi * wi;
Chris@42 111 }
Chris@42 112
Chris@42 113 /* convolution: IFFT by FFT with real/imag input/output swapped */
Chris@42 114 {
Chris@42 115 plan_dft *cldf = (plan_dft *)ego->cldf;
Chris@42 116 cldf->apply(ego->cldf, b, b+1, b, b+1);
Chris@42 117 }
Chris@42 118
Chris@42 119 /* multiply output by conjugate bluestein sequence */
Chris@42 120 for (i = 0; i < n; ++i) {
Chris@42 121 E xi = b[2*i], xr = b[2*i+1];
Chris@42 122 E wr = w[2*i], wi = w[2*i+1];
Chris@42 123 ro[i*os] = xr * wr + xi * wi;
Chris@42 124 io[i*os] = xi * wr - xr * wi;
Chris@42 125 }
Chris@42 126
Chris@42 127 X(ifree)(b);
Chris@42 128 }
Chris@42 129
Chris@42 130 static void awake(plan *ego_, enum wakefulness wakefulness)
Chris@42 131 {
Chris@42 132 P *ego = (P *) ego_;
Chris@42 133
Chris@42 134 X(plan_awake)(ego->cldf, wakefulness);
Chris@42 135
Chris@42 136 switch (wakefulness) {
Chris@42 137 case SLEEPY:
Chris@42 138 X(ifree0)(ego->w); ego->w = 0;
Chris@42 139 X(ifree0)(ego->W); ego->W = 0;
Chris@42 140 break;
Chris@42 141 default:
Chris@42 142 A(!ego->w);
Chris@42 143 mktwiddle(wakefulness, ego);
Chris@42 144 break;
Chris@42 145 }
Chris@42 146 }
Chris@42 147
Chris@42 148 static int applicable(const solver *ego, const problem *p_,
Chris@42 149 const planner *plnr)
Chris@42 150 {
Chris@42 151 const problem_dft *p = (const problem_dft *) p_;
Chris@42 152 UNUSED(ego);
Chris@42 153 return (1
Chris@42 154 && p->sz->rnk == 1
Chris@42 155 && p->vecsz->rnk == 0
Chris@42 156 /* FIXME: allow other sizes */
Chris@42 157 && X(is_prime)(p->sz->dims[0].n)
Chris@42 158
Chris@42 159 /* FIXME: avoid infinite recursion of bluestein with itself.
Chris@42 160 This works because all factors in child problems are 2, 3, 5 */
Chris@42 161 && p->sz->dims[0].n > 16
Chris@42 162
Chris@42 163 && CIMPLIES(NO_SLOWP(plnr), p->sz->dims[0].n > BLUESTEIN_MAX_SLOW)
Chris@42 164 );
Chris@42 165 }
Chris@42 166
Chris@42 167 static void destroy(plan *ego_)
Chris@42 168 {
Chris@42 169 P *ego = (P *) ego_;
Chris@42 170 X(plan_destroy_internal)(ego->cldf);
Chris@42 171 }
Chris@42 172
Chris@42 173 static void print(const plan *ego_, printer *p)
Chris@42 174 {
Chris@42 175 const P *ego = (const P *)ego_;
Chris@42 176 p->print(p, "(dft-bluestein-%D/%D%(%p%))",
Chris@42 177 ego->n, ego->nb, ego->cldf);
Chris@42 178 }
Chris@42 179
Chris@42 180 static INT choose_transform_size(INT minsz)
Chris@42 181 {
Chris@42 182 while (!X(factors_into_small_primes)(minsz))
Chris@42 183 ++minsz;
Chris@42 184 return minsz;
Chris@42 185 }
Chris@42 186
Chris@42 187 static plan *mkplan(const solver *ego, const problem *p_, planner *plnr)
Chris@42 188 {
Chris@42 189 const problem_dft *p = (const problem_dft *) p_;
Chris@42 190 P *pln;
Chris@42 191 INT n, nb;
Chris@42 192 plan *cldf = 0;
Chris@42 193 R *buf = (R *) 0;
Chris@42 194
Chris@42 195 static const plan_adt padt = {
Chris@42 196 X(dft_solve), awake, print, destroy
Chris@42 197 };
Chris@42 198
Chris@42 199 if (!applicable(ego, p_, plnr))
Chris@42 200 return (plan *) 0;
Chris@42 201
Chris@42 202 n = p->sz->dims[0].n;
Chris@42 203 nb = choose_transform_size(2 * n - 1);
Chris@42 204 buf = (R *) MALLOC(2 * nb * sizeof(R), BUFFERS);
Chris@42 205
Chris@42 206 cldf = X(mkplan_f_d)(plnr,
Chris@42 207 X(mkproblem_dft_d)(X(mktensor_1d)(nb, 2, 2),
Chris@42 208 X(mktensor_1d)(1, 0, 0),
Chris@42 209 buf, buf+1,
Chris@42 210 buf, buf+1),
Chris@42 211 NO_SLOW, 0, 0);
Chris@42 212 if (!cldf) goto nada;
Chris@42 213
Chris@42 214 X(ifree)(buf);
Chris@42 215
Chris@42 216 pln = MKPLAN_DFT(P, &padt, apply);
Chris@42 217
Chris@42 218 pln->n = n;
Chris@42 219 pln->nb = nb;
Chris@42 220 pln->w = 0;
Chris@42 221 pln->W = 0;
Chris@42 222 pln->cldf = cldf;
Chris@42 223 pln->is = p->sz->dims[0].is;
Chris@42 224 pln->os = p->sz->dims[0].os;
Chris@42 225
Chris@42 226 X(ops_add)(&cldf->ops, &cldf->ops, &pln->super.super.ops);
Chris@42 227 pln->super.super.ops.add += 4 * n + 2 * nb;
Chris@42 228 pln->super.super.ops.mul += 8 * n + 4 * nb;
Chris@42 229 pln->super.super.ops.other += 6 * (n + nb);
Chris@42 230
Chris@42 231 return &(pln->super.super);
Chris@42 232
Chris@42 233 nada:
Chris@42 234 X(ifree0)(buf);
Chris@42 235 X(plan_destroy_internal)(cldf);
Chris@42 236 return (plan *)0;
Chris@42 237 }
Chris@42 238
Chris@42 239
Chris@42 240 static solver *mksolver(void)
Chris@42 241 {
Chris@42 242 static const solver_adt sadt = { PROBLEM_DFT, mkplan, 0 };
Chris@42 243 S *slv = MKSOLVER(S, &sadt);
Chris@42 244 return &(slv->super);
Chris@42 245 }
Chris@42 246
Chris@42 247 void X(dft_bluestein_register)(planner *p)
Chris@42 248 {
Chris@42 249 REGISTER_SOLVER(p, mksolver());
Chris@42 250 }