Chris@42: /*
Chris@42:  * Copyright (c) 2003, 2007-14 Matteo Frigo
Chris@42:  * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
Chris@42:  *
Chris@42:  * This program is free software; you can redistribute it and/or modify
Chris@42:  * it under the terms of the GNU General Public License as published by
Chris@42:  * the Free Software Foundation; either version 2 of the License, or
Chris@42:  * (at your option) any later version.
Chris@42:  *
Chris@42:  * This program is distributed in the hope that it will be useful,
Chris@42:  * but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@42:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Chris@42:  * GNU General Public License for more details.
Chris@42:  *
Chris@42:  * You should have received a copy of the GNU General Public License
Chris@42:  * along with this program; if not, write to the Free Software
Chris@42:  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
Chris@42:  *
Chris@42:  */
Chris@42: 
Chris@42: 
Chris@42: #include "ct-hc2c.h"
Chris@42: 
Chris@42: typedef struct {
Chris@42:      hc2c_solver super;
Chris@42:      const hc2c_desc *desc;
Chris@42:      int bufferedp;
Chris@42:      khc2c k;
Chris@42: } S;
Chris@42: 
Chris@42: typedef struct {
Chris@42:      plan_hc2c super;
Chris@42:      khc2c k;
Chris@42:      plan *cld0, *cldm; /* children for 0th and middle butterflies */
Chris@42:      INT r, m, v, extra_iter;
Chris@42:      INT ms, vs;
Chris@42:      stride rs, brs;
Chris@42:      twid *td;
Chris@42:      const S *slv;
Chris@42: } P;
Chris@42: 
Chris@42: /*************************************************************
Chris@42:   Nonbuffered code
Chris@42:  *************************************************************/
Chris@42: static void apply(const plan *ego_, R *cr, R *ci)
Chris@42: {
Chris@42:      const P *ego = (const P *) ego_;
Chris@42:      plan_rdft2 *cld0 = (plan_rdft2 *) ego->cld0;
Chris@42:      plan_rdft2 *cldm = (plan_rdft2 *) ego->cldm;
Chris@42:      INT i, m = ego->m, v = ego->v;
Chris@42:      INT ms = ego->ms, vs = ego->vs;
Chris@42: 
Chris@42:      for (i = 0; i < v; ++i, cr += vs, ci += vs) {
Chris@42: 	  cld0->apply((plan *) cld0, cr, ci, cr, ci);
Chris@42: 	  ego->k(cr + ms, ci + ms, cr + (m-1)*ms, ci + (m-1)*ms,
Chris@42: 		 ego->td->W, ego->rs, 1, (m+1)/2, ms);
Chris@42: 	  cldm->apply((plan *) cldm, cr + (m/2)*ms, ci + (m/2)*ms, 
Chris@42: 		      cr + (m/2)*ms, ci + (m/2)*ms);
Chris@42:      }
Chris@42: }
Chris@42: 
Chris@42: static void apply_extra_iter(const plan *ego_, R *cr, R *ci)
Chris@42: {
Chris@42:      const P *ego = (const P *) ego_;
Chris@42:      plan_rdft2 *cld0 = (plan_rdft2 *) ego->cld0;
Chris@42:      plan_rdft2 *cldm = (plan_rdft2 *) ego->cldm;
Chris@42:      INT i, m = ego->m, v = ego->v;
Chris@42:      INT ms = ego->ms, vs = ego->vs;
Chris@42:      INT mm = (m-1)/2;
Chris@42: 
Chris@42:      for (i = 0; i < v; ++i, cr += vs, ci += vs) {
Chris@42: 	  cld0->apply((plan *) cld0, cr, ci, cr, ci);
Chris@42: 
Chris@42: 	  /* for 4-way SIMD when (m+1)/2-1 is odd: iterate over an
Chris@42: 	     even vector length MM-1, and then execute the last
Chris@42: 	     iteration as a 2-vector with vector stride 0.  The
Chris@42: 	     twiddle factors of the second half of the last iteration
Chris@42: 	     are bogus, but we only store the results of the first
Chris@42: 	     half. */
Chris@42: 	  ego->k(cr + ms, ci + ms, cr + (m-1)*ms, ci + (m-1)*ms,
Chris@42: 		 ego->td->W, ego->rs, 1, mm, ms);
Chris@42: 	  ego->k(cr + mm*ms, ci + mm*ms, cr + (m-mm)*ms, ci + (m-mm)*ms,
Chris@42: 		 ego->td->W, ego->rs, mm, mm+2, 0);
Chris@42: 	  cldm->apply((plan *) cldm, cr + (m/2)*ms, ci + (m/2)*ms, 
Chris@42: 		      cr + (m/2)*ms, ci + (m/2)*ms);
Chris@42:      }
Chris@42: 
Chris@42: }
Chris@42: 
Chris@42: /*************************************************************
Chris@42:   Buffered code
Chris@42:  *************************************************************/
Chris@42: 
Chris@42: /* should not be 2^k to avoid associativity conflicts */
Chris@42: static INT compute_batchsize(INT radix)
Chris@42: {
Chris@42:      /* round up to multiple of 4 */
Chris@42:      radix += 3;
Chris@42:      radix &= -4;
Chris@42: 
Chris@42:      return (radix + 2);
Chris@42: }
Chris@42: 
Chris@42: static void dobatch(const P *ego, R *Rp, R *Ip, R *Rm, R *Im,
Chris@42: 		    INT mb, INT me, INT extra_iter, R *bufp)
Chris@42: {
Chris@42:      INT b = WS(ego->brs, 1);
Chris@42:      INT rs = WS(ego->rs, 1);
Chris@42:      INT ms = ego->ms;
Chris@42:      R *bufm = bufp + b - 2;
Chris@42:      INT n = me - mb;
Chris@42: 
Chris@42:      X(cpy2d_pair_ci)(Rp + mb * ms, Ip + mb * ms, bufp, bufp + 1,
Chris@42: 		      ego->r / 2, rs, b,
Chris@42: 		      n, ms, 2);
Chris@42:      X(cpy2d_pair_ci)(Rm - mb * ms, Im - mb * ms, bufm, bufm + 1,
Chris@42: 		      ego->r / 2, rs, b,
Chris@42: 		      n, -ms, -2);
Chris@42: 
Chris@42:      if (extra_iter) {
Chris@42:           /* initialize the extra_iter element to 0.  It would be ok
Chris@42:              to leave it uninitialized, since we transform uninitialized
Chris@42:              data and ignore the result.  However, we want to avoid
Chris@42:              FP exceptions in case somebody is trapping them. */
Chris@42:           A(n < compute_batchsize(ego->r));
Chris@42:           X(zero1d_pair)(bufp + 2*n, bufp + 1 + 2*n, ego->r / 2, b);
Chris@42:           X(zero1d_pair)(bufm - 2*n, bufm + 1 - 2*n, ego->r / 2, b);
Chris@42:      }
Chris@42: 
Chris@42:      ego->k(bufp, bufp + 1, bufm, bufm + 1, ego->td->W, 
Chris@42: 	    ego->brs, mb, me + extra_iter, 2);
Chris@42:      X(cpy2d_pair_co)(bufp, bufp + 1, Rp + mb * ms, Ip + mb * ms, 
Chris@42: 		      ego->r / 2, b, rs,
Chris@42: 		      n, 2, ms);
Chris@42:      X(cpy2d_pair_co)(bufm, bufm + 1, Rm - mb * ms, Im - mb * ms,
Chris@42: 		      ego->r / 2, b, rs,
Chris@42: 		      n, -2, -ms);
Chris@42: }
Chris@42: 
Chris@42: static void apply_buf(const plan *ego_, R *cr, R *ci)
Chris@42: {
Chris@42:      const P *ego = (const P *) ego_;
Chris@42:      plan_rdft2 *cld0 = (plan_rdft2 *) ego->cld0;
Chris@42:      plan_rdft2 *cldm = (plan_rdft2 *) ego->cldm;
Chris@42:      INT i, j, ms = ego->ms, v = ego->v;
Chris@42:      INT batchsz = compute_batchsize(ego->r);
Chris@42:      R *buf;
Chris@42:      INT mb = 1, me = (ego->m+1) / 2;
Chris@42:      size_t bufsz = ego->r * batchsz * 2 * sizeof(R);
Chris@42: 
Chris@42:      BUF_ALLOC(R *, buf, bufsz);
Chris@42: 
Chris@42:      for (i = 0; i < v; ++i, cr += ego->vs, ci += ego->vs) {
Chris@42: 	  R *Rp = cr;
Chris@42: 	  R *Ip = ci;
Chris@42: 	  R *Rm = cr + ego->m * ms;
Chris@42: 	  R *Im = ci + ego->m * ms;
Chris@42: 
Chris@42: 	  cld0->apply((plan *) cld0, Rp, Ip, Rp, Ip);
Chris@42: 
Chris@42: 	  for (j = mb; j + batchsz < me; j += batchsz) 
Chris@42: 	       dobatch(ego, Rp, Ip, Rm, Im, j, j + batchsz, 0, buf);
Chris@42: 
Chris@42: 	  dobatch(ego, Rp, Ip, Rm, Im, j, me, ego->extra_iter, buf);
Chris@42: 
Chris@42: 	  cldm->apply((plan *) cldm, 
Chris@42: 		      Rp + me * ms, Ip + me * ms,
Chris@42: 		      Rp + me * ms, Ip + me * ms);
Chris@42: 
Chris@42:      }
Chris@42: 
Chris@42:      BUF_FREE(buf, bufsz);
Chris@42: }
Chris@42: 
Chris@42: /*************************************************************
Chris@42:   common code
Chris@42:  *************************************************************/
Chris@42: static void awake(plan *ego_, enum wakefulness wakefulness)
Chris@42: {
Chris@42:      P *ego = (P *) ego_;
Chris@42: 
Chris@42:      X(plan_awake)(ego->cld0, wakefulness);
Chris@42:      X(plan_awake)(ego->cldm, wakefulness);
Chris@42:      X(twiddle_awake)(wakefulness, &ego->td, ego->slv->desc->tw, 
Chris@42: 		      ego->r * ego->m, ego->r, 
Chris@42: 		      (ego->m - 1) / 2 + ego->extra_iter);
Chris@42: }
Chris@42: 
Chris@42: static void destroy(plan *ego_)
Chris@42: {
Chris@42:      P *ego = (P *) ego_;
Chris@42:      X(plan_destroy_internal)(ego->cld0);
Chris@42:      X(plan_destroy_internal)(ego->cldm);
Chris@42:      X(stride_destroy)(ego->rs);
Chris@42:      X(stride_destroy)(ego->brs);
Chris@42: }
Chris@42: 
Chris@42: static void print(const plan *ego_, printer *p)
Chris@42: {
Chris@42:      const P *ego = (const P *) ego_;
Chris@42:      const S *slv = ego->slv;
Chris@42:      const hc2c_desc *e = slv->desc;
Chris@42: 
Chris@42:      if (slv->bufferedp)
Chris@42: 	  p->print(p, "(hc2c-directbuf/%D-%D/%D/%D%v \"%s\"%(%p%)%(%p%))",
Chris@42: 		   compute_batchsize(ego->r),
Chris@42: 		   ego->r, X(twiddle_length)(ego->r, e->tw),
Chris@42: 		   ego->extra_iter, ego->v, e->nam, 
Chris@42: 		   ego->cld0, ego->cldm);
Chris@42:      else
Chris@42: 	  p->print(p, "(hc2c-direct-%D/%D/%D%v \"%s\"%(%p%)%(%p%))",
Chris@42: 		   ego->r, X(twiddle_length)(ego->r, e->tw), 
Chris@42: 		   ego->extra_iter, ego->v, e->nam, 
Chris@42: 		   ego->cld0, ego->cldm);
Chris@42: }
Chris@42: 
Chris@42: static int applicable0(const S *ego, rdft_kind kind,
Chris@42: 		       INT r, INT rs,
Chris@42: 		       INT m, INT ms, 
Chris@42: 		       INT v, INT vs,
Chris@42: 		       const R *cr, const R *ci,
Chris@42: 		       const planner *plnr,
Chris@42: 		       INT *extra_iter)
Chris@42: {
Chris@42:      const hc2c_desc *e = ego->desc;
Chris@42:      UNUSED(v);
Chris@42: 
Chris@42:      return (
Chris@42: 	  1
Chris@42: 	  && r == e->radix
Chris@42: 	  && kind == e->genus->kind
Chris@42: 
Chris@42: 	  /* first v-loop iteration */
Chris@42: 	  && ((*extra_iter = 0,
Chris@42: 	       e->genus->okp(cr + ms, ci + ms, cr + (m-1)*ms, ci + (m-1)*ms,
Chris@42: 			     rs, 1, (m+1)/2, ms, plnr))
Chris@42:               ||
Chris@42: 	      (*extra_iter = 1,
Chris@42: 	       ((e->genus->okp(cr + ms, ci + ms, cr + (m-1)*ms, ci + (m-1)*ms,
Chris@42: 			       rs, 1, (m-1)/2, ms, plnr))
Chris@42: 		&&
Chris@42: 		(e->genus->okp(cr + ms, ci + ms, cr + (m-1)*ms, ci + (m-1)*ms,
Chris@42: 			       rs, (m-1)/2, (m-1)/2 + 2, 0, plnr)))))
Chris@42: 	  
Chris@42: 	  /* subsequent v-loop iterations */
Chris@42: 	  && (cr += vs, ci += vs, 1)
Chris@42: 
Chris@42: 	  && e->genus->okp(cr + ms, ci + ms, cr + (m-1)*ms, ci + (m-1)*ms,
Chris@42: 			   rs, 1, (m+1)/2 - *extra_iter, ms, plnr)
Chris@42: 	  );
Chris@42: }
Chris@42: 
Chris@42: static int applicable0_buf(const S *ego, rdft_kind kind,
Chris@42: 			   INT r, INT rs,
Chris@42: 			   INT m, INT ms, 
Chris@42: 			   INT v, INT vs,
Chris@42: 			   const R *cr, const R *ci,
Chris@42: 			   const planner *plnr, INT *extra_iter)
Chris@42: {
Chris@42:      const hc2c_desc *e = ego->desc;
Chris@42:      INT batchsz, brs;
Chris@42:      UNUSED(v); UNUSED(rs); UNUSED(ms); UNUSED(vs);
Chris@42: 
Chris@42:      return (
Chris@42: 	  1
Chris@42: 	  && r == e->radix
Chris@42: 	  && kind == e->genus->kind
Chris@42: 
Chris@42: 	  /* ignore cr, ci, use buffer */
Chris@42: 	  && (cr = (const R *)0, ci = cr + 1, 
Chris@42: 	      batchsz = compute_batchsize(r), 
Chris@42: 	      brs = 4 * batchsz, 1)
Chris@42: 
Chris@42: 	  && e->genus->okp(cr, ci, cr + brs - 2, ci + brs - 2, 
Chris@42: 			   brs, 1, 1+batchsz, 2, plnr)
Chris@42: 
Chris@42: 	  && ((*extra_iter = 0,
Chris@42: 	       e->genus->okp(cr, ci, cr + brs - 2, ci + brs - 2, 
Chris@42: 			     brs, 1, 1 + (((m-1)/2) % batchsz), 2, plnr))
Chris@42: 	      ||
Chris@42: 	      (*extra_iter = 1,
Chris@42: 	       e->genus->okp(cr, ci, cr + brs - 2, ci + brs - 2, 
Chris@42: 			     brs, 1, 1 + 1 + (((m-1)/2) % batchsz), 2, plnr)))
Chris@42: 	      
Chris@42: 	  );
Chris@42: }
Chris@42: 
Chris@42: static int applicable(const S *ego, rdft_kind kind,
Chris@42: 		      INT r, INT rs,
Chris@42: 		      INT m, INT ms, 
Chris@42: 		      INT v, INT vs,
Chris@42: 		      R *cr, R *ci,
Chris@42: 		      const planner *plnr, INT *extra_iter)
Chris@42: {
Chris@42:      if (ego->bufferedp) {
Chris@42: 	  if (!applicable0_buf(ego, kind, r, rs, m, ms, v, vs, cr, ci, plnr,
Chris@42: 			       extra_iter))
Chris@42: 	       return 0;
Chris@42:      } else {
Chris@42: 	  if (!applicable0(ego, kind, r, rs, m, ms, v, vs, cr, ci, plnr,
Chris@42: 			   extra_iter))
Chris@42: 	       return 0;
Chris@42:      }
Chris@42: 
Chris@42:      if (NO_UGLYP(plnr) && X(ct_uglyp)((ego->bufferedp? (INT)512 : (INT)16),
Chris@42: 				       v, m * r, r))
Chris@42: 	  return 0;
Chris@42: 
Chris@42:      return 1;
Chris@42: }
Chris@42: 
Chris@42: static plan *mkcldw(const hc2c_solver *ego_, rdft_kind kind,
Chris@42: 		    INT r, INT rs,
Chris@42: 		    INT m, INT ms, 
Chris@42: 		    INT v, INT vs,
Chris@42: 		    R *cr, R *ci,
Chris@42: 		    planner *plnr)
Chris@42: {
Chris@42:      const S *ego = (const S *) ego_;
Chris@42:      P *pln;
Chris@42:      const hc2c_desc *e = ego->desc;
Chris@42:      plan *cld0 = 0, *cldm = 0;
Chris@42:      INT imid = (m / 2) * ms;
Chris@42:      INT extra_iter;
Chris@42: 
Chris@42:      static const plan_adt padt = {
Chris@42: 	  0, awake, print, destroy
Chris@42:      };
Chris@42: 
Chris@42:      if (!applicable(ego, kind, r, rs, m, ms, v, vs, cr, ci, plnr, 
Chris@42: 		     &extra_iter))
Chris@42:           return (plan *)0;
Chris@42: 
Chris@42:      cld0 = X(mkplan_d)(
Chris@42: 	  plnr, 
Chris@42: 	  X(mkproblem_rdft2_d)(X(mktensor_1d)(r, rs, rs),
Chris@42: 			       X(mktensor_0d)(),
Chris@42: 			       TAINT(cr, vs), TAINT(ci, vs),
Chris@42: 			       TAINT(cr, vs), TAINT(ci, vs),
Chris@42: 			       kind));
Chris@42:      if (!cld0) goto nada;
Chris@42: 
Chris@42:      cldm = X(mkplan_d)(
Chris@42: 	  plnr, 
Chris@42: 	  X(mkproblem_rdft2_d)(((m % 2) ?
Chris@42: 				X(mktensor_0d)() : X(mktensor_1d)(r, rs, rs) ),
Chris@42: 			       X(mktensor_0d)(),
Chris@42: 			       TAINT(cr + imid, vs), TAINT(ci + imid, vs),
Chris@42: 			       TAINT(cr + imid, vs), TAINT(ci + imid, vs),
Chris@42: 			       kind == R2HC ? R2HCII : HC2RIII));
Chris@42:      if (!cldm) goto nada;
Chris@42: 
Chris@42:      if (ego->bufferedp)
Chris@42: 	  pln = MKPLAN_HC2C(P, &padt, apply_buf);
Chris@42:      else
Chris@42: 	  pln = MKPLAN_HC2C(P, &padt, extra_iter ? apply_extra_iter : apply);
Chris@42: 
Chris@42:      pln->k = ego->k;
Chris@42:      pln->td = 0;
Chris@42:      pln->r = r; pln->rs = X(mkstride)(r, rs);
Chris@42:      pln->m = m; pln->ms = ms;
Chris@42:      pln->v = v; pln->vs = vs;
Chris@42:      pln->slv = ego;
Chris@42:      pln->brs = X(mkstride)(r, 4 * compute_batchsize(r));
Chris@42:      pln->cld0 = cld0;
Chris@42:      pln->cldm = cldm;
Chris@42:      pln->extra_iter = extra_iter;
Chris@42: 
Chris@42:      X(ops_zero)(&pln->super.super.ops);
Chris@42:      X(ops_madd2)(v * (((m - 1) / 2) / e->genus->vl),
Chris@42: 		  &e->ops, &pln->super.super.ops);
Chris@42:      X(ops_madd2)(v, &cld0->ops, &pln->super.super.ops);
Chris@42:      X(ops_madd2)(v, &cldm->ops, &pln->super.super.ops);
Chris@42: 
Chris@42:      if (ego->bufferedp) 
Chris@42: 	  pln->super.super.ops.other += 4 * r * m * v;
Chris@42: 
Chris@42:      return &(pln->super.super);
Chris@42: 
Chris@42:  nada:
Chris@42:      X(plan_destroy_internal)(cld0);
Chris@42:      X(plan_destroy_internal)(cldm);
Chris@42:      return 0;
Chris@42: }
Chris@42: 
Chris@42: static void regone(planner *plnr, khc2c codelet,
Chris@42: 		   const hc2c_desc *desc, 
Chris@42: 		   hc2c_kind hc2ckind, 
Chris@42: 		   int bufferedp)
Chris@42: {
Chris@42:      S *slv = (S *)X(mksolver_hc2c)(sizeof(S), desc->radix, hc2ckind, mkcldw);
Chris@42:      slv->k = codelet;
Chris@42:      slv->desc = desc;
Chris@42:      slv->bufferedp = bufferedp;
Chris@42:      REGISTER_SOLVER(plnr, &(slv->super.super));
Chris@42: }
Chris@42: 
Chris@42: void X(regsolver_hc2c_direct)(planner *plnr, khc2c codelet,
Chris@42: 			      const hc2c_desc *desc,
Chris@42: 			      hc2c_kind hc2ckind)
Chris@42: {
Chris@42:      regone(plnr, codelet, desc, hc2ckind, /* bufferedp */0);
Chris@42:      regone(plnr, codelet, desc, hc2ckind, /* bufferedp */1);
Chris@42: }