Mercurial > hg > sv-dependency-builds

annotate src/fftw-3.3.5/dft/ct.c @ 42:2cd0e3b3e1fd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Current fftw source
author Chris Cannam
date Tue, 18 Oct 2016 13:40:26 +0100
parents
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
Chris@42 22 #include "ct.h"
Chris@42 23
Chris@42 24 ct_solver *(*X(mksolver_ct_hook))(size_t, INT, int,
Chris@42 25 ct_mkinferior, ct_force_vrecursion) = 0;
Chris@42 26
Chris@42 27 typedef struct {
Chris@42 28 plan_dft super;
Chris@42 29 plan *cld;
Chris@42 30 plan *cldw;
Chris@42 31 INT r;
Chris@42 32 } P;
Chris@42 33
Chris@42 34 static void apply_dit(const plan *ego_, R *ri, R *ii, R *ro, R *io)
Chris@42 35 {
Chris@42 36 const P *ego = (const P *) ego_;
Chris@42 37 plan_dft *cld;
Chris@42 38 plan_dftw *cldw;
Chris@42 39
Chris@42 40 cld = (plan_dft *) ego->cld;
Chris@42 41 cld->apply(ego->cld, ri, ii, ro, io);
Chris@42 42
Chris@42 43 cldw = (plan_dftw *) ego->cldw;
Chris@42 44 cldw->apply(ego->cldw, ro, io);
Chris@42 45 }
Chris@42 46
Chris@42 47 static void apply_dif(const plan *ego_, R *ri, R *ii, R *ro, R *io)
Chris@42 48 {
Chris@42 49 const P *ego = (const P *) ego_;
Chris@42 50 plan_dft *cld;
Chris@42 51 plan_dftw *cldw;
Chris@42 52
Chris@42 53 cldw = (plan_dftw *) ego->cldw;
Chris@42 54 cldw->apply(ego->cldw, ri, ii);
Chris@42 55
Chris@42 56 cld = (plan_dft *) ego->cld;
Chris@42 57 cld->apply(ego->cld, ri, ii, ro, io);
Chris@42 58 }
Chris@42 59
Chris@42 60 static void awake(plan *ego_, enum wakefulness wakefulness)
Chris@42 61 {
Chris@42 62 P *ego = (P *) ego_;
Chris@42 63 X(plan_awake)(ego->cld, wakefulness);
Chris@42 64 X(plan_awake)(ego->cldw, wakefulness);
Chris@42 65 }
Chris@42 66
Chris@42 67 static void destroy(plan *ego_)
Chris@42 68 {
Chris@42 69 P *ego = (P *) ego_;
Chris@42 70 X(plan_destroy_internal)(ego->cldw);
Chris@42 71 X(plan_destroy_internal)(ego->cld);
Chris@42 72 }
Chris@42 73
Chris@42 74 static void print(const plan *ego_, printer *p)
Chris@42 75 {
Chris@42 76 const P *ego = (const P *) ego_;
Chris@42 77 p->print(p, "(dft-ct-%s/%D%(%p%)%(%p%))",
Chris@42 78 ego->super.apply == apply_dit ? "dit" : "dif",
Chris@42 79 ego->r, ego->cldw, ego->cld);
Chris@42 80 }
Chris@42 81
Chris@42 82 static int applicable0(const ct_solver *ego, const problem *p_, planner *plnr)
Chris@42 83 {
Chris@42 84 const problem_dft *p = (const problem_dft *) p_;
Chris@42 85 INT r;
Chris@42 86
Chris@42 87 return (1
Chris@42 88 && p->sz->rnk == 1
Chris@42 89 && p->vecsz->rnk <= 1
Chris@42 90
Chris@42 91 /* DIF destroys the input and we don't like it */
Chris@42 92 && (ego->dec == DECDIT ||
Chris@42 93 p->ri == p->ro ||
Chris@42 94 !NO_DESTROY_INPUTP(plnr))
Chris@42 95
Chris@42 96 && ((r = X(choose_radix)(ego->r, p->sz->dims[0].n)) > 1)
Chris@42 97 && p->sz->dims[0].n > r);
Chris@42 98 }
Chris@42 99
Chris@42 100
Chris@42 101 int X(ct_applicable)(const ct_solver *ego, const problem *p_, planner *plnr)
Chris@42 102 {
Chris@42 103 const problem_dft *p;
Chris@42 104
Chris@42 105 if (!applicable0(ego, p_, plnr))
Chris@42 106 return 0;
Chris@42 107
Chris@42 108 p = (const problem_dft *) p_;
Chris@42 109
Chris@42 110 return (0
Chris@42 111 || ego->dec == DECDIF+TRANSPOSE
Chris@42 112 || p->vecsz->rnk == 0
Chris@42 113 || !NO_VRECURSEP(plnr)
Chris@42 114 || (ego->force_vrecursionp && ego->force_vrecursionp(ego, p))
Chris@42 115 );
Chris@42 116 }
Chris@42 117
Chris@42 118
Chris@42 119 static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
Chris@42 120 {
Chris@42 121 const ct_solver *ego = (const ct_solver *) ego_;
Chris@42 122 const problem_dft *p;
Chris@42 123 P *pln = 0;
Chris@42 124 plan *cld = 0, *cldw = 0;
Chris@42 125 INT n, r, m, v, ivs, ovs;
Chris@42 126 iodim *d;
Chris@42 127
Chris@42 128 static const plan_adt padt = {
Chris@42 129 X(dft_solve), awake, print, destroy
Chris@42 130 };
Chris@42 131
Chris@42 132 if ((NO_NONTHREADEDP(plnr)) || !X(ct_applicable)(ego, p_, plnr))
Chris@42 133 return (plan *) 0;
Chris@42 134
Chris@42 135 p = (const problem_dft *) p_;
Chris@42 136 d = p->sz->dims;
Chris@42 137 n = d[0].n;
Chris@42 138 r = X(choose_radix)(ego->r, n);
Chris@42 139 m = n / r;
Chris@42 140
Chris@42 141 X(tensor_tornk1)(p->vecsz, &v, &ivs, &ovs);
Chris@42 142
Chris@42 143 switch (ego->dec) {
Chris@42 144 case DECDIT:
Chris@42 145 {
Chris@42 146 cldw = ego->mkcldw(ego,
Chris@42 147 r, m * d[0].os, m * d[0].os,
Chris@42 148 m, d[0].os,
Chris@42 149 v, ovs, ovs,
Chris@42 150 0, m,
Chris@42 151 p->ro, p->io, plnr);
Chris@42 152 if (!cldw) goto nada;
Chris@42 153
Chris@42 154 cld = X(mkplan_d)(plnr,
Chris@42 155 X(mkproblem_dft_d)(
Chris@42 156 X(mktensor_1d)(m, r * d[0].is, d[0].os),
Chris@42 157 X(mktensor_2d)(r, d[0].is, m * d[0].os,
Chris@42 158 v, ivs, ovs),
Chris@42 159 p->ri, p->ii, p->ro, p->io)
Chris@42 160 );
Chris@42 161 if (!cld) goto nada;
Chris@42 162
Chris@42 163 pln = MKPLAN_DFT(P, &padt, apply_dit);
Chris@42 164 break;
Chris@42 165 }
Chris@42 166 case DECDIF:
Chris@42 167 case DECDIF+TRANSPOSE:
Chris@42 168 {
Chris@42 169 INT cors, covs; /* cldw ors, ovs */
Chris@42 170 if (ego->dec == DECDIF+TRANSPOSE) {
Chris@42 171 cors = ivs;
Chris@42 172 covs = m * d[0].is;
Chris@42 173 /* ensure that we generate well-formed dftw subproblems */
Chris@42 174 /* FIXME: too conservative */
Chris@42 175 if (!(1
Chris@42 176 && r == v
Chris@42 177 && d[0].is == r * cors))
Chris@42 178 goto nada;
Chris@42 179
Chris@42 180 /* FIXME: allow in-place only for now, like in
Chris@42 181 fftw-3.[01] */
Chris@42 182 if (!(1
Chris@42 183 && p->ri == p->ro
Chris@42 184 && d[0].is == r * d[0].os
Chris@42 185 && cors == d[0].os
Chris@42 186 && covs == ovs
Chris@42 187 ))
Chris@42 188 goto nada;
Chris@42 189 } else {
Chris@42 190 cors = m * d[0].is;
Chris@42 191 covs = ivs;
Chris@42 192 }
Chris@42 193
Chris@42 194 cldw = ego->mkcldw(ego,
Chris@42 195 r, m * d[0].is, cors,
Chris@42 196 m, d[0].is,
Chris@42 197 v, ivs, covs,
Chris@42 198 0, m,
Chris@42 199 p->ri, p->ii, plnr);
Chris@42 200 if (!cldw) goto nada;
Chris@42 201
Chris@42 202 cld = X(mkplan_d)(plnr,
Chris@42 203 X(mkproblem_dft_d)(
Chris@42 204 X(mktensor_1d)(m, d[0].is, r * d[0].os),
Chris@42 205 X(mktensor_2d)(r, cors, d[0].os,
Chris@42 206 v, covs, ovs),
Chris@42 207 p->ri, p->ii, p->ro, p->io)
Chris@42 208 );
Chris@42 209 if (!cld) goto nada;
Chris@42 210
Chris@42 211 pln = MKPLAN_DFT(P, &padt, apply_dif);
Chris@42 212 break;
Chris@42 213 }
Chris@42 214
Chris@42 215 default: A(0);
Chris@42 216
Chris@42 217 }
Chris@42 218
Chris@42 219 pln->cld = cld;
Chris@42 220 pln->cldw = cldw;
Chris@42 221 pln->r = r;
Chris@42 222 X(ops_add)(&cld->ops, &cldw->ops, &pln->super.super.ops);
Chris@42 223
Chris@42 224 /* inherit could_prune_now_p attribute from cldw */
Chris@42 225 pln->super.super.could_prune_now_p = cldw->could_prune_now_p;
Chris@42 226 return &(pln->super.super);
Chris@42 227
Chris@42 228 nada:
Chris@42 229 X(plan_destroy_internal)(cldw);
Chris@42 230 X(plan_destroy_internal)(cld);
Chris@42 231 return (plan *) 0;
Chris@42 232 }
Chris@42 233
Chris@42 234 ct_solver *X(mksolver_ct)(size_t size, INT r, int dec,
Chris@42 235 ct_mkinferior mkcldw,
Chris@42 236 ct_force_vrecursion force_vrecursionp)
Chris@42 237 {
Chris@42 238 static const solver_adt sadt = { PROBLEM_DFT, mkplan, 0 };
Chris@42 239 ct_solver *slv = (ct_solver *)X(mksolver)(size, &sadt);
Chris@42 240 slv->r = r;
Chris@42 241 slv->dec = dec;
Chris@42 242 slv->mkcldw = mkcldw;
Chris@42 243 slv->force_vrecursionp = force_vrecursionp;
Chris@42 244 return slv;
Chris@42 245 }
Chris@42 246
Chris@42 247 plan *X(mkplan_dftw)(size_t size, const plan_adt *adt, dftwapply apply)
Chris@42 248 {
Chris@42 249 plan_dftw *ego;
Chris@42 250
Chris@42 251 ego = (plan_dftw *) X(mkplan)(size, adt);
Chris@42 252 ego->apply = apply;
Chris@42 253
Chris@42 254 return &(ego->super);
Chris@42 255 }