Mercurial > hg > sv-dependency-builds

annotate src/fftw-3.3.5/rdft/ct-hc2c.c @ 82:d0c2a83c1364

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add FFTW 3.3.8 source, and a Linux build
author Chris Cannam
date Tue, 19 Nov 2019 14:52:55 +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 "ct-hc2c.h"
Chris@42 22 #include "dft.h"
Chris@42 23
Chris@42 24 typedef struct {
Chris@42 25 plan_rdft2 super;
Chris@42 26 plan *cld;
Chris@42 27 plan *cldw;
Chris@42 28 INT r;
Chris@42 29 } P;
Chris@42 30
Chris@42 31 static void apply_dit(const plan *ego_, R *r0, R *r1, R *cr, R *ci)
Chris@42 32 {
Chris@42 33 const P *ego = (const P *) ego_;
Chris@42 34 plan_rdft *cld;
Chris@42 35 plan_hc2c *cldw;
Chris@42 36 UNUSED(r1);
Chris@42 37
Chris@42 38 cld = (plan_rdft *) ego->cld;
Chris@42 39 cld->apply(ego->cld, r0, cr);
Chris@42 40
Chris@42 41 cldw = (plan_hc2c *) ego->cldw;
Chris@42 42 cldw->apply(ego->cldw, cr, ci);
Chris@42 43 }
Chris@42 44
Chris@42 45 static void apply_dif(const plan *ego_, R *r0, R *r1, R *cr, R *ci)
Chris@42 46 {
Chris@42 47 const P *ego = (const P *) ego_;
Chris@42 48 plan_rdft *cld;
Chris@42 49 plan_hc2c *cldw;
Chris@42 50 UNUSED(r1);
Chris@42 51
Chris@42 52 cldw = (plan_hc2c *) ego->cldw;
Chris@42 53 cldw->apply(ego->cldw, cr, ci);
Chris@42 54
Chris@42 55 cld = (plan_rdft *) ego->cld;
Chris@42 56 cld->apply(ego->cld, cr, r0);
Chris@42 57 }
Chris@42 58
Chris@42 59 static void apply_dit_dft(const plan *ego_, R *r0, R *r1, R *cr, R *ci)
Chris@42 60 {
Chris@42 61 const P *ego = (const P *) ego_;
Chris@42 62 plan_dft *cld;
Chris@42 63 plan_hc2c *cldw;
Chris@42 64
Chris@42 65 cld = (plan_dft *) ego->cld;
Chris@42 66 cld->apply(ego->cld, r0, r1, cr, ci);
Chris@42 67
Chris@42 68 cldw = (plan_hc2c *) ego->cldw;
Chris@42 69 cldw->apply(ego->cldw, cr, ci);
Chris@42 70 }
Chris@42 71
Chris@42 72 static void apply_dif_dft(const plan *ego_, R *r0, R *r1, R *cr, R *ci)
Chris@42 73 {
Chris@42 74 const P *ego = (const P *) ego_;
Chris@42 75 plan_dft *cld;
Chris@42 76 plan_hc2c *cldw;
Chris@42 77
Chris@42 78 cldw = (plan_hc2c *) ego->cldw;
Chris@42 79 cldw->apply(ego->cldw, cr, ci);
Chris@42 80
Chris@42 81 cld = (plan_dft *) ego->cld;
Chris@42 82 cld->apply(ego->cld, ci, cr, r1, r0);
Chris@42 83 }
Chris@42 84
Chris@42 85 static void awake(plan *ego_, enum wakefulness wakefulness)
Chris@42 86 {
Chris@42 87 P *ego = (P *) ego_;
Chris@42 88 X(plan_awake)(ego->cld, wakefulness);
Chris@42 89 X(plan_awake)(ego->cldw, wakefulness);
Chris@42 90 }
Chris@42 91
Chris@42 92 static void destroy(plan *ego_)
Chris@42 93 {
Chris@42 94 P *ego = (P *) ego_;
Chris@42 95 X(plan_destroy_internal)(ego->cldw);
Chris@42 96 X(plan_destroy_internal)(ego->cld);
Chris@42 97 }
Chris@42 98
Chris@42 99 static void print(const plan *ego_, printer *p)
Chris@42 100 {
Chris@42 101 const P *ego = (const P *) ego_;
Chris@42 102 p->print(p, "(rdft2-ct-%s/%D%(%p%)%(%p%))",
Chris@42 103 (ego->super.apply == apply_dit ||
Chris@42 104 ego->super.apply == apply_dit_dft)
Chris@42 105 ? "dit" : "dif",
Chris@42 106 ego->r, ego->cldw, ego->cld);
Chris@42 107 }
Chris@42 108
Chris@42 109 static int applicable0(const hc2c_solver *ego, const problem *p_, planner *plnr)
Chris@42 110 {
Chris@42 111 const problem_rdft2 *p = (const problem_rdft2 *) p_;
Chris@42 112 INT r;
Chris@42 113
Chris@42 114 return (1
Chris@42 115 && p->sz->rnk == 1
Chris@42 116 && p->vecsz->rnk <= 1
Chris@42 117
Chris@42 118 && (/* either the problem is R2HC, which is solved by DIT */
Chris@42 119 (p->kind == R2HC)
Chris@42 120 ||
Chris@42 121 /* or the problem is HC2R, in which case it is solved
Chris@42 122 by DIF, which destroys the input */
Chris@42 123 (p->kind == HC2R &&
Chris@42 124 (p->r0 == p->cr || !NO_DESTROY_INPUTP(plnr))))
Chris@42 125
Chris@42 126 && ((r = X(choose_radix)(ego->r, p->sz->dims[0].n)) > 0)
Chris@42 127 && p->sz->dims[0].n > r);
Chris@42 128 }
Chris@42 129
Chris@42 130 static int hc2c_applicable(const hc2c_solver *ego, const problem *p_,
Chris@42 131 planner *plnr)
Chris@42 132 {
Chris@42 133 const problem_rdft2 *p;
Chris@42 134
Chris@42 135 if (!applicable0(ego, p_, plnr))
Chris@42 136 return 0;
Chris@42 137
Chris@42 138 p = (const problem_rdft2 *) p_;
Chris@42 139
Chris@42 140 return (0
Chris@42 141 || p->vecsz->rnk == 0
Chris@42 142 || !NO_VRECURSEP(plnr)
Chris@42 143 );
Chris@42 144 }
Chris@42 145
Chris@42 146 static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
Chris@42 147 {
Chris@42 148 const hc2c_solver *ego = (const hc2c_solver *) ego_;
Chris@42 149 const problem_rdft2 *p;
Chris@42 150 P *pln = 0;
Chris@42 151 plan *cld = 0, *cldw = 0;
Chris@42 152 INT n, r, m, v, ivs, ovs;
Chris@42 153 iodim *d;
Chris@42 154
Chris@42 155 static const plan_adt padt = {
Chris@42 156 X(rdft2_solve), awake, print, destroy
Chris@42 157 };
Chris@42 158
Chris@42 159 if (!hc2c_applicable(ego, p_, plnr))
Chris@42 160 return (plan *) 0;
Chris@42 161
Chris@42 162 p = (const problem_rdft2 *) p_;
Chris@42 163 d = p->sz->dims;
Chris@42 164 n = d[0].n;
Chris@42 165 r = X(choose_radix)(ego->r, n);
Chris@42 166 A((r % 2) == 0);
Chris@42 167 m = n / r;
Chris@42 168
Chris@42 169 X(tensor_tornk1)(p->vecsz, &v, &ivs, &ovs);
Chris@42 170
Chris@42 171 switch (p->kind) {
Chris@42 172 case R2HC:
Chris@42 173 cldw = ego->mkcldw(ego, R2HC,
Chris@42 174 r, m * d[0].os,
Chris@42 175 m, d[0].os,
Chris@42 176 v, ovs,
Chris@42 177 p->cr, p->ci, plnr);
Chris@42 178 if (!cldw) goto nada;
Chris@42 179
Chris@42 180 switch (ego->hc2ckind) {
Chris@42 181 case HC2C_VIA_RDFT:
Chris@42 182 cld = X(mkplan_d)(
Chris@42 183 plnr,
Chris@42 184 X(mkproblem_rdft_1_d)(
Chris@42 185 X(mktensor_1d)(m, (r/2)*d[0].is, d[0].os),
Chris@42 186 X(mktensor_3d)(
Chris@42 187 2, p->r1 - p->r0, p->ci - p->cr,
Chris@42 188 r / 2, d[0].is, m * d[0].os,
Chris@42 189 v, ivs, ovs),
Chris@42 190 p->r0, p->cr, R2HC)
Chris@42 191 );
Chris@42 192 if (!cld) goto nada;
Chris@42 193
Chris@42 194 pln = MKPLAN_RDFT2(P, &padt, apply_dit);
Chris@42 195 break;
Chris@42 196
Chris@42 197 case HC2C_VIA_DFT:
Chris@42 198 cld = X(mkplan_d)(
Chris@42 199 plnr,
Chris@42 200 X(mkproblem_dft_d)(
Chris@42 201 X(mktensor_1d)(m, (r/2)*d[0].is, d[0].os),
Chris@42 202 X(mktensor_2d)(
Chris@42 203 r / 2, d[0].is, m * d[0].os,
Chris@42 204 v, ivs, ovs),
Chris@42 205 p->r0, p->r1, p->cr, p->ci)
Chris@42 206 );
Chris@42 207 if (!cld) goto nada;
Chris@42 208
Chris@42 209 pln = MKPLAN_RDFT2(P, &padt, apply_dit_dft);
Chris@42 210 break;
Chris@42 211 }
Chris@42 212 break;
Chris@42 213
Chris@42 214 case HC2R:
Chris@42 215 cldw = ego->mkcldw(ego, HC2R,
Chris@42 216 r, m * d[0].is,
Chris@42 217 m, d[0].is,
Chris@42 218 v, ivs,
Chris@42 219 p->cr, p->ci, plnr);
Chris@42 220 if (!cldw) goto nada;
Chris@42 221
Chris@42 222 switch (ego->hc2ckind) {
Chris@42 223 case HC2C_VIA_RDFT:
Chris@42 224 cld = X(mkplan_d)(
Chris@42 225 plnr,
Chris@42 226 X(mkproblem_rdft_1_d)(
Chris@42 227 X(mktensor_1d)(m, d[0].is, (r/2)*d[0].os),
Chris@42 228 X(mktensor_3d)(
Chris@42 229 2, p->ci - p->cr, p->r1 - p->r0,
Chris@42 230 r / 2, m * d[0].is, d[0].os,
Chris@42 231 v, ivs, ovs),
Chris@42 232 p->cr, p->r0, HC2R)
Chris@42 233 );
Chris@42 234 if (!cld) goto nada;
Chris@42 235
Chris@42 236 pln = MKPLAN_RDFT2(P, &padt, apply_dif);
Chris@42 237 break;
Chris@42 238
Chris@42 239 case HC2C_VIA_DFT:
Chris@42 240 cld = X(mkplan_d)(
Chris@42 241 plnr,
Chris@42 242 X(mkproblem_dft_d)(
Chris@42 243 X(mktensor_1d)(m, d[0].is, (r/2)*d[0].os),
Chris@42 244 X(mktensor_2d)(
Chris@42 245 r / 2, m * d[0].is, d[0].os,
Chris@42 246 v, ivs, ovs),
Chris@42 247 p->ci, p->cr, p->r1, p->r0)
Chris@42 248 );
Chris@42 249 if (!cld) goto nada;
Chris@42 250
Chris@42 251 pln = MKPLAN_RDFT2(P, &padt, apply_dif_dft);
Chris@42 252 break;
Chris@42 253 }
Chris@42 254 break;
Chris@42 255
Chris@42 256 default:
Chris@42 257 A(0);
Chris@42 258 }
Chris@42 259
Chris@42 260 pln->cld = cld;
Chris@42 261 pln->cldw = cldw;
Chris@42 262 pln->r = r;
Chris@42 263 X(ops_add)(&cld->ops, &cldw->ops, &pln->super.super.ops);
Chris@42 264
Chris@42 265 /* inherit could_prune_now_p attribute from cldw */
Chris@42 266 pln->super.super.could_prune_now_p = cldw->could_prune_now_p;
Chris@42 267
Chris@42 268 return &(pln->super.super);
Chris@42 269
Chris@42 270 nada:
Chris@42 271 X(plan_destroy_internal)(cldw);
Chris@42 272 X(plan_destroy_internal)(cld);
Chris@42 273 return (plan *) 0;
Chris@42 274 }
Chris@42 275
Chris@42 276 hc2c_solver *X(mksolver_hc2c)(size_t size, INT r,
Chris@42 277 hc2c_kind hc2ckind,
Chris@42 278 hc2c_mkinferior mkcldw)
Chris@42 279 {
Chris@42 280 static const solver_adt sadt = { PROBLEM_RDFT2, mkplan, 0 };
Chris@42 281 hc2c_solver *slv = (hc2c_solver *)X(mksolver)(size, &sadt);
Chris@42 282 slv->r = r;
Chris@42 283 slv->hc2ckind = hc2ckind;
Chris@42 284 slv->mkcldw = mkcldw;
Chris@42 285 return slv;
Chris@42 286 }
Chris@42 287
Chris@42 288 plan *X(mkplan_hc2c)(size_t size, const plan_adt *adt, hc2capply apply)
Chris@42 289 {
Chris@42 290 plan_hc2c *ego;
Chris@42 291
Chris@42 292 ego = (plan_hc2c *) X(mkplan)(size, adt);
Chris@42 293 ego->apply = apply;
Chris@42 294
Chris@42 295 return &(ego->super);
Chris@42 296 }