Chris@10: /* Chris@10: * Copyright (c) 2003, 2007-11 Matteo Frigo Chris@10: * Copyright (c) 2003, 2007-11 Massachusetts Institute of Technology Chris@10: * Chris@10: * This program is free software; you can redistribute it and/or modify Chris@10: * it under the terms of the GNU General Public License as published by Chris@10: * the Free Software Foundation; either version 2 of the License, or Chris@10: * (at your option) any later version. Chris@10: * Chris@10: * This program is distributed in the hope that it will be useful, Chris@10: * but WITHOUT ANY WARRANTY; without even the implied warranty of Chris@10: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the Chris@10: * GNU General Public License for more details. Chris@10: * Chris@10: * You should have received a copy of the GNU General Public License Chris@10: * along with this program; if not, write to the Free Software Chris@10: * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA Chris@10: * Chris@10: */ Chris@10: Chris@10: Chris@10: #include "hc2hc.h" Chris@10: Chris@10: typedef struct { Chris@10: hc2hc_solver super; Chris@10: const hc2hc_desc *desc; Chris@10: khc2hc k; Chris@10: int bufferedp; Chris@10: } S; Chris@10: Chris@10: typedef struct { Chris@10: plan_hc2hc super; Chris@10: khc2hc k; Chris@10: plan *cld0, *cldm; /* children for 0th and middle butterflies */ Chris@10: INT r, m, v; Chris@10: INT ms, vs, mb, me; Chris@10: stride rs, brs; Chris@10: twid *td; Chris@10: const S *slv; Chris@10: } P; Chris@10: Chris@10: /************************************************************* Chris@10: Nonbuffered code Chris@10: *************************************************************/ Chris@10: static void apply(const plan *ego_, R *IO) Chris@10: { Chris@10: const P *ego = (const P *) ego_; Chris@10: plan_rdft *cld0 = (plan_rdft *) ego->cld0; Chris@10: plan_rdft *cldm = (plan_rdft *) ego->cldm; Chris@10: INT i, m = ego->m, v = ego->v; Chris@10: INT mb = ego->mb, me = ego->me; Chris@10: INT ms = ego->ms, vs = ego->vs; Chris@10: Chris@10: for (i = 0; i < v; ++i, IO += vs) { Chris@10: cld0->apply((plan *) cld0, IO, IO); Chris@10: ego->k(IO + ms * mb, IO + (m - mb) * ms, Chris@10: ego->td->W, ego->rs, mb, me, ms); Chris@10: cldm->apply((plan *) cldm, IO + (m/2) * ms, IO + (m/2) * ms); Chris@10: } Chris@10: } Chris@10: Chris@10: /************************************************************* Chris@10: Buffered code Chris@10: *************************************************************/ Chris@10: Chris@10: /* should not be 2^k to avoid associativity conflicts */ Chris@10: static INT compute_batchsize(INT radix) Chris@10: { Chris@10: /* round up to multiple of 4 */ Chris@10: radix += 3; Chris@10: radix &= -4; Chris@10: Chris@10: return (radix + 2); Chris@10: } Chris@10: Chris@10: static void dobatch(const P *ego, R *IOp, R *IOm, Chris@10: INT mb, INT me, R *bufp) Chris@10: { Chris@10: INT b = WS(ego->brs, 1); Chris@10: INT rs = WS(ego->rs, 1); Chris@10: INT r = ego->r; Chris@10: INT ms = ego->ms; Chris@10: R *bufm = bufp + b - 1; Chris@10: Chris@10: X(cpy2d_ci)(IOp + mb * ms, bufp, r, rs, b, me - mb, ms, 1, 1); Chris@10: X(cpy2d_ci)(IOm - mb * ms, bufm, r, rs, b, me - mb, -ms, -1, 1); Chris@10: Chris@10: ego->k(bufp, bufm, ego->td->W, ego->brs, mb, me, 1); Chris@10: Chris@10: X(cpy2d_co)(bufp, IOp + mb * ms, r, b, rs, me - mb, 1, ms, 1); Chris@10: X(cpy2d_co)(bufm, IOm - mb * ms, r, b, rs, me - mb, -1, -ms, 1); Chris@10: } Chris@10: Chris@10: static void apply_buf(const plan *ego_, R *IO) Chris@10: { Chris@10: const P *ego = (const P *) ego_; Chris@10: plan_rdft *cld0 = (plan_rdft *) ego->cld0; Chris@10: plan_rdft *cldm = (plan_rdft *) ego->cldm; Chris@10: INT i, j, m = ego->m, v = ego->v, r = ego->r; Chris@10: INT mb = ego->mb, me = ego->me, ms = ego->ms; Chris@10: INT batchsz = compute_batchsize(r); Chris@10: R *buf; Chris@10: size_t bufsz = r * batchsz * 2 * sizeof(R); Chris@10: Chris@10: BUF_ALLOC(R *, buf, bufsz); Chris@10: Chris@10: for (i = 0; i < v; ++i, IO += ego->vs) { Chris@10: R *IOp = IO; Chris@10: R *IOm = IO + m * ms; Chris@10: Chris@10: cld0->apply((plan *) cld0, IO, IO); Chris@10: Chris@10: for (j = mb; j + batchsz < me; j += batchsz) Chris@10: dobatch(ego, IOp, IOm, j, j + batchsz, buf); Chris@10: Chris@10: dobatch(ego, IOp, IOm, j, me, buf); Chris@10: Chris@10: cldm->apply((plan *) cldm, IO + ms * (m/2), IO + ms * (m/2)); Chris@10: } Chris@10: Chris@10: BUF_FREE(buf, bufsz); Chris@10: } Chris@10: Chris@10: static void awake(plan *ego_, enum wakefulness wakefulness) Chris@10: { Chris@10: P *ego = (P *) ego_; Chris@10: Chris@10: X(plan_awake)(ego->cld0, wakefulness); Chris@10: X(plan_awake)(ego->cldm, wakefulness); Chris@10: X(twiddle_awake)(wakefulness, &ego->td, ego->slv->desc->tw, Chris@10: ego->r * ego->m, ego->r, (ego->m - 1) / 2); Chris@10: } Chris@10: Chris@10: static void destroy(plan *ego_) Chris@10: { Chris@10: P *ego = (P *) ego_; Chris@10: X(plan_destroy_internal)(ego->cld0); Chris@10: X(plan_destroy_internal)(ego->cldm); Chris@10: X(stride_destroy)(ego->rs); Chris@10: X(stride_destroy)(ego->brs); Chris@10: } Chris@10: Chris@10: static void print(const plan *ego_, printer *p) Chris@10: { Chris@10: const P *ego = (const P *) ego_; Chris@10: const S *slv = ego->slv; Chris@10: const hc2hc_desc *e = slv->desc; Chris@10: INT batchsz = compute_batchsize(ego->r); Chris@10: Chris@10: if (slv->bufferedp) Chris@10: p->print(p, "(hc2hc-directbuf/%D-%D/%D%v \"%s\"%(%p%)%(%p%))", Chris@10: batchsz, ego->r, X(twiddle_length)(ego->r, e->tw), Chris@10: ego->v, e->nam, ego->cld0, ego->cldm); Chris@10: else Chris@10: p->print(p, "(hc2hc-direct-%D/%D%v \"%s\"%(%p%)%(%p%))", Chris@10: ego->r, X(twiddle_length)(ego->r, e->tw), ego->v, e->nam, Chris@10: ego->cld0, ego->cldm); Chris@10: } Chris@10: Chris@10: static int applicable0(const S *ego, rdft_kind kind, INT r) Chris@10: { Chris@10: const hc2hc_desc *e = ego->desc; Chris@10: Chris@10: return (1 Chris@10: && r == e->radix Chris@10: && kind == e->genus->kind Chris@10: ); Chris@10: } Chris@10: Chris@10: static int applicable(const S *ego, rdft_kind kind, INT r, INT m, INT v, Chris@10: const planner *plnr) Chris@10: { Chris@10: if (!applicable0(ego, kind, r)) Chris@10: return 0; Chris@10: Chris@10: if (NO_UGLYP(plnr) && X(ct_uglyp)((ego->bufferedp? (INT)512 : (INT)16), Chris@10: v, m * r, r)) Chris@10: return 0; Chris@10: Chris@10: return 1; Chris@10: } Chris@10: Chris@10: #define CLDMP(m, mstart, mcount) (2 * ((mstart) + (mcount)) == (m) + 2) Chris@10: #define CLD0P(mstart) ((mstart) == 0) Chris@10: Chris@10: static plan *mkcldw(const hc2hc_solver *ego_, Chris@10: rdft_kind kind, INT r, INT m, INT ms, INT v, INT vs, Chris@10: INT mstart, INT mcount, Chris@10: R *IO, planner *plnr) Chris@10: { Chris@10: const S *ego = (const S *) ego_; Chris@10: P *pln; Chris@10: const hc2hc_desc *e = ego->desc; Chris@10: plan *cld0 = 0, *cldm = 0; Chris@10: INT imid = (m / 2) * ms; Chris@10: INT rs = m * ms; Chris@10: Chris@10: static const plan_adt padt = { Chris@10: 0, awake, print, destroy Chris@10: }; Chris@10: Chris@10: if (!applicable(ego, kind, r, m, v, plnr)) Chris@10: return (plan *)0; Chris@10: Chris@10: cld0 = X(mkplan_d)( Chris@10: plnr, Chris@10: X(mkproblem_rdft_1_d)((CLD0P(mstart) ? Chris@10: X(mktensor_1d)(r, rs, rs) : X(mktensor_0d)()), Chris@10: X(mktensor_0d)(), Chris@10: TAINT(IO, vs), TAINT(IO, vs), Chris@10: kind)); Chris@10: if (!cld0) goto nada; Chris@10: Chris@10: cldm = X(mkplan_d)( Chris@10: plnr, Chris@10: X(mkproblem_rdft_1_d)((CLDMP(m, mstart, mcount) ? Chris@10: X(mktensor_1d)(r, rs, rs) : X(mktensor_0d)()), Chris@10: X(mktensor_0d)(), Chris@10: TAINT(IO + imid, vs), TAINT(IO + imid, vs), Chris@10: kind == R2HC ? R2HCII : HC2RIII)); Chris@10: if (!cldm) goto nada; Chris@10: Chris@10: pln = MKPLAN_HC2HC(P, &padt, ego->bufferedp ? apply_buf : apply); Chris@10: Chris@10: pln->k = ego->k; Chris@10: pln->td = 0; Chris@10: pln->r = r; pln->rs = X(mkstride)(r, rs); Chris@10: pln->m = m; pln->ms = ms; Chris@10: pln->v = v; pln->vs = vs; Chris@10: pln->slv = ego; Chris@10: pln->brs = X(mkstride)(r, 2 * compute_batchsize(r)); Chris@10: pln->cld0 = cld0; Chris@10: pln->cldm = cldm; Chris@10: pln->mb = mstart + CLD0P(mstart); Chris@10: pln->me = mstart + mcount - CLDMP(m, mstart, mcount); Chris@10: Chris@10: X(ops_zero)(&pln->super.super.ops); Chris@10: X(ops_madd2)(v * ((pln->me - pln->mb) / e->genus->vl), Chris@10: &e->ops, &pln->super.super.ops); Chris@10: X(ops_madd2)(v, &cld0->ops, &pln->super.super.ops); Chris@10: X(ops_madd2)(v, &cldm->ops, &pln->super.super.ops); Chris@10: Chris@10: if (ego->bufferedp) Chris@10: pln->super.super.ops.other += 4 * r * (pln->me - pln->mb) * v; Chris@10: Chris@10: pln->super.super.could_prune_now_p = Chris@10: (!ego->bufferedp && r >= 5 && r < 64 && m >= r); Chris@10: Chris@10: return &(pln->super.super); Chris@10: Chris@10: nada: Chris@10: X(plan_destroy_internal)(cld0); Chris@10: X(plan_destroy_internal)(cldm); Chris@10: return 0; Chris@10: } Chris@10: Chris@10: static void regone(planner *plnr, khc2hc codelet, const hc2hc_desc *desc, Chris@10: int bufferedp) Chris@10: { Chris@10: S *slv = (S *)X(mksolver_hc2hc)(sizeof(S), desc->radix, mkcldw); Chris@10: slv->k = codelet; Chris@10: slv->desc = desc; Chris@10: slv->bufferedp = bufferedp; Chris@10: REGISTER_SOLVER(plnr, &(slv->super.super)); Chris@10: if (X(mksolver_hc2hc_hook)) { Chris@10: slv = (S *)X(mksolver_hc2hc_hook)(sizeof(S), desc->radix, mkcldw); Chris@10: slv->k = codelet; Chris@10: slv->desc = desc; Chris@10: slv->bufferedp = bufferedp; Chris@10: REGISTER_SOLVER(plnr, &(slv->super.super)); Chris@10: } Chris@10: } Chris@10: Chris@10: void X(regsolver_hc2hc_direct)(planner *plnr, khc2hc codelet, Chris@10: const hc2hc_desc *desc) Chris@10: { Chris@10: regone(plnr, codelet, desc, /* bufferedp */0); Chris@10: regone(plnr, codelet, desc, /* bufferedp */1); Chris@10: }