Mercurial > hg > sv-dependency-builds

annotate src/fftw-3.3.5/dft/bluestein.c @ 83:ae30d91d2ffe

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Replace these with versions built using an older toolset (so as to avoid ABI compatibilities when linking on Ubuntu 14.04 for packaging purposes)
author Chris Cannam
date Fri, 07 Feb 2020 11:51:13 +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 }