cannam@127: /* cannam@127: * Copyright (c) 2003, 2007-14 Matteo Frigo cannam@127: * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology cannam@127: * cannam@127: * This program is free software; you can redistribute it and/or modify cannam@127: * it under the terms of the GNU General Public License as published by cannam@127: * the Free Software Foundation; either version 2 of the License, or cannam@127: * (at your option) any later version. cannam@127: * cannam@127: * This program is distributed in the hope that it will be useful, cannam@127: * but WITHOUT ANY WARRANTY; without even the implied warranty of cannam@127: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the cannam@127: * GNU General Public License for more details. cannam@127: * cannam@127: * You should have received a copy of the GNU General Public License cannam@127: * along with this program; if not, write to the Free Software cannam@127: * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA cannam@127: * cannam@127: */ cannam@127: cannam@127: /* express a twiddle problem in terms of dft + multiplication by cannam@127: twiddle factors */ cannam@127: cannam@127: #include "ct.h" cannam@127: cannam@127: typedef ct_solver S; cannam@127: cannam@127: typedef struct { cannam@127: plan_dftw super; cannam@127: cannam@127: INT r, rs, m, mb, me, ms, v, vs; cannam@127: cannam@127: plan *cld; cannam@127: cannam@127: twid *td; cannam@127: cannam@127: const S *slv; cannam@127: int dec; cannam@127: } P; cannam@127: cannam@127: static void mktwiddle(P *ego, enum wakefulness wakefulness) cannam@127: { cannam@127: static const tw_instr tw[] = { { TW_FULL, 0, 0 }, { TW_NEXT, 1, 0 } }; cannam@127: cannam@127: /* note that R and M are swapped, to allow for sequential cannam@127: access both to data and twiddles */ cannam@127: X(twiddle_awake)(wakefulness, &ego->td, tw, cannam@127: ego->r * ego->m, ego->m, ego->r); cannam@127: } cannam@127: cannam@127: static void bytwiddle(const P *ego, R *rio, R *iio) cannam@127: { cannam@127: INT iv, ir, im; cannam@127: INT r = ego->r, rs = ego->rs; cannam@127: INT m = ego->m, mb = ego->mb, me = ego->me, ms = ego->ms; cannam@127: INT v = ego->v, vs = ego->vs; cannam@127: const R *W = ego->td->W; cannam@127: cannam@127: mb += (mb == 0); /* skip m=0 iteration */ cannam@127: for (iv = 0; iv < v; ++iv) { cannam@127: for (ir = 1; ir < r; ++ir) { cannam@127: for (im = mb; im < me; ++im) { cannam@127: R *pr = rio + ms * im + rs * ir; cannam@127: R *pi = iio + ms * im + rs * ir; cannam@127: E xr = *pr; cannam@127: E xi = *pi; cannam@127: E wr = W[2 * im + (2 * (m-1)) * ir - 2]; cannam@127: E wi = W[2 * im + (2 * (m-1)) * ir - 1]; cannam@127: *pr = xr * wr + xi * wi; cannam@127: *pi = xi * wr - xr * wi; cannam@127: } cannam@127: } cannam@127: rio += vs; cannam@127: iio += vs; cannam@127: } cannam@127: } cannam@127: cannam@127: static int applicable(INT irs, INT ors, INT ivs, INT ovs, cannam@127: const planner *plnr) cannam@127: { cannam@127: return (1 cannam@127: && irs == ors cannam@127: && ivs == ovs cannam@127: && !NO_SLOWP(plnr) cannam@127: ); cannam@127: } cannam@127: cannam@127: static void apply_dit(const plan *ego_, R *rio, R *iio) cannam@127: { cannam@127: const P *ego = (const P *) ego_; cannam@127: plan_dft *cld; cannam@127: INT dm = ego->ms * ego->mb; cannam@127: cannam@127: bytwiddle(ego, rio, iio); cannam@127: cannam@127: cld = (plan_dft *) ego->cld; cannam@127: cld->apply(ego->cld, rio + dm, iio + dm, rio + dm, iio + dm); cannam@127: } cannam@127: cannam@127: static void apply_dif(const plan *ego_, R *rio, R *iio) cannam@127: { cannam@127: const P *ego = (const P *) ego_; cannam@127: plan_dft *cld; cannam@127: INT dm = ego->ms * ego->mb; cannam@127: cannam@127: cld = (plan_dft *) ego->cld; cannam@127: cld->apply(ego->cld, rio + dm, iio + dm, rio + dm, iio + dm); cannam@127: cannam@127: bytwiddle(ego, rio, iio); cannam@127: } cannam@127: cannam@127: static void awake(plan *ego_, enum wakefulness wakefulness) cannam@127: { cannam@127: P *ego = (P *) ego_; cannam@127: X(plan_awake)(ego->cld, wakefulness); cannam@127: mktwiddle(ego, wakefulness); cannam@127: } cannam@127: cannam@127: static void destroy(plan *ego_) cannam@127: { cannam@127: P *ego = (P *) ego_; cannam@127: X(plan_destroy_internal)(ego->cld); cannam@127: } cannam@127: cannam@127: static void print(const plan *ego_, printer *p) cannam@127: { cannam@127: const P *ego = (const P *) ego_; cannam@127: p->print(p, "(dftw-generic-%s-%D-%D%v%(%p%))", cannam@127: ego->dec == DECDIT ? "dit" : "dif", cannam@127: ego->r, ego->m, ego->v, ego->cld); cannam@127: } cannam@127: cannam@127: static plan *mkcldw(const ct_solver *ego_, cannam@127: INT r, INT irs, INT ors, cannam@127: INT m, INT ms, cannam@127: INT v, INT ivs, INT ovs, cannam@127: INT mstart, INT mcount, cannam@127: R *rio, R *iio, cannam@127: planner *plnr) cannam@127: { cannam@127: const S *ego = (const S *)ego_; cannam@127: P *pln; cannam@127: plan *cld = 0; cannam@127: INT dm = ms * mstart; cannam@127: cannam@127: static const plan_adt padt = { cannam@127: 0, awake, print, destroy cannam@127: }; cannam@127: cannam@127: A(mstart >= 0 && mstart + mcount <= m); cannam@127: if (!applicable(irs, ors, ivs, ovs, plnr)) cannam@127: return (plan *)0; cannam@127: cannam@127: cld = X(mkplan_d)(plnr, cannam@127: X(mkproblem_dft_d)( cannam@127: X(mktensor_1d)(r, irs, irs), cannam@127: X(mktensor_2d)(mcount, ms, ms, v, ivs, ivs), cannam@127: rio + dm, iio + dm, rio + dm, iio + dm) cannam@127: ); cannam@127: if (!cld) goto nada; cannam@127: cannam@127: pln = MKPLAN_DFTW(P, &padt, ego->dec == DECDIT ? apply_dit : apply_dif); cannam@127: pln->slv = ego; cannam@127: pln->cld = cld; cannam@127: pln->r = r; cannam@127: pln->rs = irs; cannam@127: pln->m = m; cannam@127: pln->ms = ms; cannam@127: pln->v = v; cannam@127: pln->vs = ivs; cannam@127: pln->mb = mstart; cannam@127: pln->me = mstart + mcount; cannam@127: pln->dec = ego->dec; cannam@127: pln->td = 0; cannam@127: cannam@127: { cannam@127: double n0 = (r - 1) * (mcount - 1) * v; cannam@127: pln->super.super.ops = cld->ops; cannam@127: pln->super.super.ops.mul += 8 * n0; cannam@127: pln->super.super.ops.add += 4 * n0; cannam@127: pln->super.super.ops.other += 8 * n0; cannam@127: } cannam@127: return &(pln->super.super); cannam@127: cannam@127: nada: cannam@127: X(plan_destroy_internal)(cld); cannam@127: return (plan *) 0; cannam@127: } cannam@127: cannam@127: static void regsolver(planner *plnr, INT r, int dec) cannam@127: { cannam@127: S *slv = (S *)X(mksolver_ct)(sizeof(S), r, dec, mkcldw, 0); cannam@127: REGISTER_SOLVER(plnr, &(slv->super)); cannam@127: if (X(mksolver_ct_hook)) { cannam@127: slv = (S *)X(mksolver_ct_hook)(sizeof(S), r, dec, mkcldw, 0); cannam@127: REGISTER_SOLVER(plnr, &(slv->super)); cannam@127: } cannam@127: } cannam@127: cannam@127: void X(ct_generic_register)(planner *p) cannam@127: { cannam@127: regsolver(p, 0, DECDIT); cannam@127: regsolver(p, 0, DECDIF); cannam@127: }