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: #include "dft.h" cannam@127: cannam@127: typedef struct { cannam@127: solver super; cannam@127: } S; cannam@127: cannam@127: typedef struct { cannam@127: plan_dft super; cannam@127: twid *td; cannam@127: INT n, is, os; cannam@127: } P; cannam@127: cannam@127: cannam@127: static void cdot(INT n, const E *x, const R *w, cannam@127: R *or0, R *oi0, R *or1, R *oi1) cannam@127: { cannam@127: INT i; cannam@127: cannam@127: E rr = x[0], ri = 0, ir = x[1], ii = 0; cannam@127: x += 2; cannam@127: for (i = 1; i + i < n; ++i) { cannam@127: rr += x[0] * w[0]; cannam@127: ir += x[1] * w[0]; cannam@127: ri += x[2] * w[1]; cannam@127: ii += x[3] * w[1]; cannam@127: x += 4; w += 2; cannam@127: } cannam@127: *or0 = rr + ii; cannam@127: *oi0 = ir - ri; cannam@127: *or1 = rr - ii; cannam@127: *oi1 = ir + ri; cannam@127: } cannam@127: cannam@127: static void hartley(INT n, const R *xr, const R *xi, INT xs, E *o, cannam@127: R *pr, R *pi) cannam@127: { cannam@127: INT i; cannam@127: E sr, si; cannam@127: o[0] = sr = xr[0]; o[1] = si = xi[0]; o += 2; cannam@127: for (i = 1; i + i < n; ++i) { cannam@127: sr += (o[0] = xr[i * xs] + xr[(n - i) * xs]); cannam@127: si += (o[1] = xi[i * xs] + xi[(n - i) * xs]); cannam@127: o[2] = xr[i * xs] - xr[(n - i) * xs]; cannam@127: o[3] = xi[i * xs] - xi[(n - i) * xs]; cannam@127: o += 4; cannam@127: } cannam@127: *pr = sr; cannam@127: *pi = si; cannam@127: } cannam@127: cannam@127: static void apply(const plan *ego_, R *ri, R *ii, R *ro, R *io) cannam@127: { cannam@127: const P *ego = (const P *) ego_; cannam@127: INT i; cannam@127: INT n = ego->n, is = ego->is, os = ego->os; cannam@127: const R *W = ego->td->W; cannam@127: E *buf; cannam@127: size_t bufsz = n * 2 * sizeof(E); cannam@127: cannam@127: BUF_ALLOC(E *, buf, bufsz); cannam@127: hartley(n, ri, ii, is, buf, ro, io); cannam@127: cannam@127: for (i = 1; i + i < n; ++i) { cannam@127: cdot(n, buf, W, cannam@127: ro + i * os, io + i * os, cannam@127: ro + (n - i) * os, io + (n - i) * os); cannam@127: W += n - 1; cannam@127: } cannam@127: cannam@127: BUF_FREE(buf, bufsz); cannam@127: } cannam@127: cannam@127: static void awake(plan *ego_, enum wakefulness wakefulness) cannam@127: { cannam@127: P *ego = (P *) ego_; cannam@127: static const tw_instr half_tw[] = { cannam@127: { TW_HALF, 1, 0 }, cannam@127: { TW_NEXT, 1, 0 } cannam@127: }; cannam@127: cannam@127: X(twiddle_awake)(wakefulness, &ego->td, half_tw, ego->n, ego->n, cannam@127: (ego->n - 1) / 2); 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: cannam@127: p->print(p, "(dft-generic-%D)", ego->n); cannam@127: } cannam@127: cannam@127: static int applicable(const solver *ego, const problem *p_, cannam@127: const planner *plnr) cannam@127: { cannam@127: const problem_dft *p = (const problem_dft *) p_; cannam@127: UNUSED(ego); cannam@127: cannam@127: return (1 cannam@127: && p->sz->rnk == 1 cannam@127: && p->vecsz->rnk == 0 cannam@127: && (p->sz->dims[0].n % 2) == 1 cannam@127: && CIMPLIES(NO_LARGE_GENERICP(plnr), p->sz->dims[0].n < GENERIC_MIN_BAD) cannam@127: && CIMPLIES(NO_SLOWP(plnr), p->sz->dims[0].n > GENERIC_MAX_SLOW) cannam@127: && X(is_prime)(p->sz->dims[0].n) cannam@127: ); cannam@127: } cannam@127: cannam@127: static plan *mkplan(const solver *ego, const problem *p_, planner *plnr) cannam@127: { cannam@127: const problem_dft *p; cannam@127: P *pln; cannam@127: INT n; cannam@127: cannam@127: static const plan_adt padt = { cannam@127: X(dft_solve), awake, print, X(plan_null_destroy) cannam@127: }; cannam@127: cannam@127: if (!applicable(ego, p_, plnr)) cannam@127: return (plan *)0; cannam@127: cannam@127: pln = MKPLAN_DFT(P, &padt, apply); cannam@127: cannam@127: p = (const problem_dft *) p_; cannam@127: pln->n = n = p->sz->dims[0].n; cannam@127: pln->is = p->sz->dims[0].is; cannam@127: pln->os = p->sz->dims[0].os; cannam@127: pln->td = 0; cannam@127: cannam@127: pln->super.super.ops.add = (n-1) * 5; cannam@127: pln->super.super.ops.mul = 0; cannam@127: pln->super.super.ops.fma = (n-1) * (n-1) ; cannam@127: #if 0 /* these are nice pipelined sequential loads and should cost nothing */ cannam@127: pln->super.super.ops.other = (n-1)*(4 + 1 + 2 * (n-1)); /* approximate */ cannam@127: #endif cannam@127: cannam@127: return &(pln->super.super); cannam@127: } cannam@127: cannam@127: static solver *mksolver(void) cannam@127: { cannam@127: static const solver_adt sadt = { PROBLEM_DFT, mkplan, 0 }; cannam@127: S *slv = MKSOLVER(S, &sadt); cannam@127: return &(slv->super); cannam@127: } cannam@127: cannam@127: void X(dft_generic_register)(planner *p) cannam@127: { cannam@127: REGISTER_SOLVER(p, mksolver()); cannam@127: }