Mercurial > hg > sv-dependency-builds

annotate src/fftw-3.3.5/dft/ct.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
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 }