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: 
cannam@127: #include "ct-hc2c.h"
cannam@127: 
cannam@127: typedef struct {
cannam@127:      hc2c_solver super;
cannam@127:      const hc2c_desc *desc;
cannam@127:      int bufferedp;
cannam@127:      khc2c k;
cannam@127: } S;
cannam@127: 
cannam@127: typedef struct {
cannam@127:      plan_hc2c super;
cannam@127:      khc2c k;
cannam@127:      plan *cld0, *cldm; /* children for 0th and middle butterflies */
cannam@127:      INT r, m, v, extra_iter;
cannam@127:      INT ms, vs;
cannam@127:      stride rs, brs;
cannam@127:      twid *td;
cannam@127:      const S *slv;
cannam@127: } P;
cannam@127: 
cannam@127: /*************************************************************
cannam@127:   Nonbuffered code
cannam@127:  *************************************************************/
cannam@127: static void apply(const plan *ego_, R *cr, R *ci)
cannam@127: {
cannam@127:      const P *ego = (const P *) ego_;
cannam@127:      plan_rdft2 *cld0 = (plan_rdft2 *) ego->cld0;
cannam@127:      plan_rdft2 *cldm = (plan_rdft2 *) ego->cldm;
cannam@127:      INT i, m = ego->m, v = ego->v;
cannam@127:      INT ms = ego->ms, vs = ego->vs;
cannam@127: 
cannam@127:      for (i = 0; i < v; ++i, cr += vs, ci += vs) {
cannam@127: 	  cld0->apply((plan *) cld0, cr, ci, cr, ci);
cannam@127: 	  ego->k(cr + ms, ci + ms, cr + (m-1)*ms, ci + (m-1)*ms,
cannam@127: 		 ego->td->W, ego->rs, 1, (m+1)/2, ms);
cannam@127: 	  cldm->apply((plan *) cldm, cr + (m/2)*ms, ci + (m/2)*ms, 
cannam@127: 		      cr + (m/2)*ms, ci + (m/2)*ms);
cannam@127:      }
cannam@127: }
cannam@127: 
cannam@127: static void apply_extra_iter(const plan *ego_, R *cr, R *ci)
cannam@127: {
cannam@127:      const P *ego = (const P *) ego_;
cannam@127:      plan_rdft2 *cld0 = (plan_rdft2 *) ego->cld0;
cannam@127:      plan_rdft2 *cldm = (plan_rdft2 *) ego->cldm;
cannam@127:      INT i, m = ego->m, v = ego->v;
cannam@127:      INT ms = ego->ms, vs = ego->vs;
cannam@127:      INT mm = (m-1)/2;
cannam@127: 
cannam@127:      for (i = 0; i < v; ++i, cr += vs, ci += vs) {
cannam@127: 	  cld0->apply((plan *) cld0, cr, ci, cr, ci);
cannam@127: 
cannam@127: 	  /* for 4-way SIMD when (m+1)/2-1 is odd: iterate over an
cannam@127: 	     even vector length MM-1, and then execute the last
cannam@127: 	     iteration as a 2-vector with vector stride 0.  The
cannam@127: 	     twiddle factors of the second half of the last iteration
cannam@127: 	     are bogus, but we only store the results of the first
cannam@127: 	     half. */
cannam@127: 	  ego->k(cr + ms, ci + ms, cr + (m-1)*ms, ci + (m-1)*ms,
cannam@127: 		 ego->td->W, ego->rs, 1, mm, ms);
cannam@127: 	  ego->k(cr + mm*ms, ci + mm*ms, cr + (m-mm)*ms, ci + (m-mm)*ms,
cannam@127: 		 ego->td->W, ego->rs, mm, mm+2, 0);
cannam@127: 	  cldm->apply((plan *) cldm, cr + (m/2)*ms, ci + (m/2)*ms, 
cannam@127: 		      cr + (m/2)*ms, ci + (m/2)*ms);
cannam@127:      }
cannam@127: 
cannam@127: }
cannam@127: 
cannam@127: /*************************************************************
cannam@127:   Buffered code
cannam@127:  *************************************************************/
cannam@127: 
cannam@127: /* should not be 2^k to avoid associativity conflicts */
cannam@127: static INT compute_batchsize(INT radix)
cannam@127: {
cannam@127:      /* round up to multiple of 4 */
cannam@127:      radix += 3;
cannam@127:      radix &= -4;
cannam@127: 
cannam@127:      return (radix + 2);
cannam@127: }
cannam@127: 
cannam@127: static void dobatch(const P *ego, R *Rp, R *Ip, R *Rm, R *Im,
cannam@127: 		    INT mb, INT me, INT extra_iter, R *bufp)
cannam@127: {
cannam@127:      INT b = WS(ego->brs, 1);
cannam@127:      INT rs = WS(ego->rs, 1);
cannam@127:      INT ms = ego->ms;
cannam@127:      R *bufm = bufp + b - 2;
cannam@127:      INT n = me - mb;
cannam@127: 
cannam@127:      X(cpy2d_pair_ci)(Rp + mb * ms, Ip + mb * ms, bufp, bufp + 1,
cannam@127: 		      ego->r / 2, rs, b,
cannam@127: 		      n, ms, 2);
cannam@127:      X(cpy2d_pair_ci)(Rm - mb * ms, Im - mb * ms, bufm, bufm + 1,
cannam@127: 		      ego->r / 2, rs, b,
cannam@127: 		      n, -ms, -2);
cannam@127: 
cannam@127:      if (extra_iter) {
cannam@127:           /* initialize the extra_iter element to 0.  It would be ok
cannam@127:              to leave it uninitialized, since we transform uninitialized
cannam@127:              data and ignore the result.  However, we want to avoid
cannam@127:              FP exceptions in case somebody is trapping them. */
cannam@127:           A(n < compute_batchsize(ego->r));
cannam@127:           X(zero1d_pair)(bufp + 2*n, bufp + 1 + 2*n, ego->r / 2, b);
cannam@127:           X(zero1d_pair)(bufm - 2*n, bufm + 1 - 2*n, ego->r / 2, b);
cannam@127:      }
cannam@127: 
cannam@127:      ego->k(bufp, bufp + 1, bufm, bufm + 1, ego->td->W, 
cannam@127: 	    ego->brs, mb, me + extra_iter, 2);
cannam@127:      X(cpy2d_pair_co)(bufp, bufp + 1, Rp + mb * ms, Ip + mb * ms, 
cannam@127: 		      ego->r / 2, b, rs,
cannam@127: 		      n, 2, ms);
cannam@127:      X(cpy2d_pair_co)(bufm, bufm + 1, Rm - mb * ms, Im - mb * ms,
cannam@127: 		      ego->r / 2, b, rs,
cannam@127: 		      n, -2, -ms);
cannam@127: }
cannam@127: 
cannam@127: static void apply_buf(const plan *ego_, R *cr, R *ci)
cannam@127: {
cannam@127:      const P *ego = (const P *) ego_;
cannam@127:      plan_rdft2 *cld0 = (plan_rdft2 *) ego->cld0;
cannam@127:      plan_rdft2 *cldm = (plan_rdft2 *) ego->cldm;
cannam@127:      INT i, j, ms = ego->ms, v = ego->v;
cannam@127:      INT batchsz = compute_batchsize(ego->r);
cannam@127:      R *buf;
cannam@127:      INT mb = 1, me = (ego->m+1) / 2;
cannam@127:      size_t bufsz = ego->r * batchsz * 2 * sizeof(R);
cannam@127: 
cannam@127:      BUF_ALLOC(R *, buf, bufsz);
cannam@127: 
cannam@127:      for (i = 0; i < v; ++i, cr += ego->vs, ci += ego->vs) {
cannam@127: 	  R *Rp = cr;
cannam@127: 	  R *Ip = ci;
cannam@127: 	  R *Rm = cr + ego->m * ms;
cannam@127: 	  R *Im = ci + ego->m * ms;
cannam@127: 
cannam@127: 	  cld0->apply((plan *) cld0, Rp, Ip, Rp, Ip);
cannam@127: 
cannam@127: 	  for (j = mb; j + batchsz < me; j += batchsz) 
cannam@127: 	       dobatch(ego, Rp, Ip, Rm, Im, j, j + batchsz, 0, buf);
cannam@127: 
cannam@127: 	  dobatch(ego, Rp, Ip, Rm, Im, j, me, ego->extra_iter, buf);
cannam@127: 
cannam@127: 	  cldm->apply((plan *) cldm, 
cannam@127: 		      Rp + me * ms, Ip + me * ms,
cannam@127: 		      Rp + me * ms, Ip + me * ms);
cannam@127: 
cannam@127:      }
cannam@127: 
cannam@127:      BUF_FREE(buf, bufsz);
cannam@127: }
cannam@127: 
cannam@127: /*************************************************************
cannam@127:   common code
cannam@127:  *************************************************************/
cannam@127: static void awake(plan *ego_, enum wakefulness wakefulness)
cannam@127: {
cannam@127:      P *ego = (P *) ego_;
cannam@127: 
cannam@127:      X(plan_awake)(ego->cld0, wakefulness);
cannam@127:      X(plan_awake)(ego->cldm, wakefulness);
cannam@127:      X(twiddle_awake)(wakefulness, &ego->td, ego->slv->desc->tw, 
cannam@127: 		      ego->r * ego->m, ego->r, 
cannam@127: 		      (ego->m - 1) / 2 + ego->extra_iter);
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->cld0);
cannam@127:      X(plan_destroy_internal)(ego->cldm);
cannam@127:      X(stride_destroy)(ego->rs);
cannam@127:      X(stride_destroy)(ego->brs);
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:      const S *slv = ego->slv;
cannam@127:      const hc2c_desc *e = slv->desc;
cannam@127: 
cannam@127:      if (slv->bufferedp)
cannam@127: 	  p->print(p, "(hc2c-directbuf/%D-%D/%D/%D%v \"%s\"%(%p%)%(%p%))",
cannam@127: 		   compute_batchsize(ego->r),
cannam@127: 		   ego->r, X(twiddle_length)(ego->r, e->tw),
cannam@127: 		   ego->extra_iter, ego->v, e->nam, 
cannam@127: 		   ego->cld0, ego->cldm);
cannam@127:      else
cannam@127: 	  p->print(p, "(hc2c-direct-%D/%D/%D%v \"%s\"%(%p%)%(%p%))",
cannam@127: 		   ego->r, X(twiddle_length)(ego->r, e->tw), 
cannam@127: 		   ego->extra_iter, ego->v, e->nam, 
cannam@127: 		   ego->cld0, ego->cldm);
cannam@127: }
cannam@127: 
cannam@127: static int applicable0(const S *ego, rdft_kind kind,
cannam@127: 		       INT r, INT rs,
cannam@127: 		       INT m, INT ms, 
cannam@127: 		       INT v, INT vs,
cannam@127: 		       const R *cr, const R *ci,
cannam@127: 		       const planner *plnr,
cannam@127: 		       INT *extra_iter)
cannam@127: {
cannam@127:      const hc2c_desc *e = ego->desc;
cannam@127:      UNUSED(v);
cannam@127: 
cannam@127:      return (
cannam@127: 	  1
cannam@127: 	  && r == e->radix
cannam@127: 	  && kind == e->genus->kind
cannam@127: 
cannam@127: 	  /* first v-loop iteration */
cannam@127: 	  && ((*extra_iter = 0,
cannam@127: 	       e->genus->okp(cr + ms, ci + ms, cr + (m-1)*ms, ci + (m-1)*ms,
cannam@127: 			     rs, 1, (m+1)/2, ms, plnr))
cannam@127:               ||
cannam@127: 	      (*extra_iter = 1,
cannam@127: 	       ((e->genus->okp(cr + ms, ci + ms, cr + (m-1)*ms, ci + (m-1)*ms,
cannam@127: 			       rs, 1, (m-1)/2, ms, plnr))
cannam@127: 		&&
cannam@127: 		(e->genus->okp(cr + ms, ci + ms, cr + (m-1)*ms, ci + (m-1)*ms,
cannam@127: 			       rs, (m-1)/2, (m-1)/2 + 2, 0, plnr)))))
cannam@127: 	  
cannam@127: 	  /* subsequent v-loop iterations */
cannam@127: 	  && (cr += vs, ci += vs, 1)
cannam@127: 
cannam@127: 	  && e->genus->okp(cr + ms, ci + ms, cr + (m-1)*ms, ci + (m-1)*ms,
cannam@127: 			   rs, 1, (m+1)/2 - *extra_iter, ms, plnr)
cannam@127: 	  );
cannam@127: }
cannam@127: 
cannam@127: static int applicable0_buf(const S *ego, rdft_kind kind,
cannam@127: 			   INT r, INT rs,
cannam@127: 			   INT m, INT ms, 
cannam@127: 			   INT v, INT vs,
cannam@127: 			   const R *cr, const R *ci,
cannam@127: 			   const planner *plnr, INT *extra_iter)
cannam@127: {
cannam@127:      const hc2c_desc *e = ego->desc;
cannam@127:      INT batchsz, brs;
cannam@127:      UNUSED(v); UNUSED(rs); UNUSED(ms); UNUSED(vs);
cannam@127: 
cannam@127:      return (
cannam@127: 	  1
cannam@127: 	  && r == e->radix
cannam@127: 	  && kind == e->genus->kind
cannam@127: 
cannam@127: 	  /* ignore cr, ci, use buffer */
cannam@127: 	  && (cr = (const R *)0, ci = cr + 1, 
cannam@127: 	      batchsz = compute_batchsize(r), 
cannam@127: 	      brs = 4 * batchsz, 1)
cannam@127: 
cannam@127: 	  && e->genus->okp(cr, ci, cr + brs - 2, ci + brs - 2, 
cannam@127: 			   brs, 1, 1+batchsz, 2, plnr)
cannam@127: 
cannam@127: 	  && ((*extra_iter = 0,
cannam@127: 	       e->genus->okp(cr, ci, cr + brs - 2, ci + brs - 2, 
cannam@127: 			     brs, 1, 1 + (((m-1)/2) % batchsz), 2, plnr))
cannam@127: 	      ||
cannam@127: 	      (*extra_iter = 1,
cannam@127: 	       e->genus->okp(cr, ci, cr + brs - 2, ci + brs - 2, 
cannam@127: 			     brs, 1, 1 + 1 + (((m-1)/2) % batchsz), 2, plnr)))
cannam@127: 	      
cannam@127: 	  );
cannam@127: }
cannam@127: 
cannam@127: static int applicable(const S *ego, rdft_kind kind,
cannam@127: 		      INT r, INT rs,
cannam@127: 		      INT m, INT ms, 
cannam@127: 		      INT v, INT vs,
cannam@127: 		      R *cr, R *ci,
cannam@127: 		      const planner *plnr, INT *extra_iter)
cannam@127: {
cannam@127:      if (ego->bufferedp) {
cannam@127: 	  if (!applicable0_buf(ego, kind, r, rs, m, ms, v, vs, cr, ci, plnr,
cannam@127: 			       extra_iter))
cannam@127: 	       return 0;
cannam@127:      } else {
cannam@127: 	  if (!applicable0(ego, kind, r, rs, m, ms, v, vs, cr, ci, plnr,
cannam@127: 			   extra_iter))
cannam@127: 	       return 0;
cannam@127:      }
cannam@127: 
cannam@127:      if (NO_UGLYP(plnr) && X(ct_uglyp)((ego->bufferedp? (INT)512 : (INT)16),
cannam@127: 				       v, m * r, r))
cannam@127: 	  return 0;
cannam@127: 
cannam@127:      return 1;
cannam@127: }
cannam@127: 
cannam@127: static plan *mkcldw(const hc2c_solver *ego_, rdft_kind kind,
cannam@127: 		    INT r, INT rs,
cannam@127: 		    INT m, INT ms, 
cannam@127: 		    INT v, INT vs,
cannam@127: 		    R *cr, R *ci,
cannam@127: 		    planner *plnr)
cannam@127: {
cannam@127:      const S *ego = (const S *) ego_;
cannam@127:      P *pln;
cannam@127:      const hc2c_desc *e = ego->desc;
cannam@127:      plan *cld0 = 0, *cldm = 0;
cannam@127:      INT imid = (m / 2) * ms;
cannam@127:      INT extra_iter;
cannam@127: 
cannam@127:      static const plan_adt padt = {
cannam@127: 	  0, awake, print, destroy
cannam@127:      };
cannam@127: 
cannam@127:      if (!applicable(ego, kind, r, rs, m, ms, v, vs, cr, ci, plnr, 
cannam@127: 		     &extra_iter))
cannam@127:           return (plan *)0;
cannam@127: 
cannam@127:      cld0 = X(mkplan_d)(
cannam@127: 	  plnr, 
cannam@127: 	  X(mkproblem_rdft2_d)(X(mktensor_1d)(r, rs, rs),
cannam@127: 			       X(mktensor_0d)(),
cannam@127: 			       TAINT(cr, vs), TAINT(ci, vs),
cannam@127: 			       TAINT(cr, vs), TAINT(ci, vs),
cannam@127: 			       kind));
cannam@127:      if (!cld0) goto nada;
cannam@127: 
cannam@127:      cldm = X(mkplan_d)(
cannam@127: 	  plnr, 
cannam@127: 	  X(mkproblem_rdft2_d)(((m % 2) ?
cannam@127: 				X(mktensor_0d)() : X(mktensor_1d)(r, rs, rs) ),
cannam@127: 			       X(mktensor_0d)(),
cannam@127: 			       TAINT(cr + imid, vs), TAINT(ci + imid, vs),
cannam@127: 			       TAINT(cr + imid, vs), TAINT(ci + imid, vs),
cannam@127: 			       kind == R2HC ? R2HCII : HC2RIII));
cannam@127:      if (!cldm) goto nada;
cannam@127: 
cannam@127:      if (ego->bufferedp)
cannam@127: 	  pln = MKPLAN_HC2C(P, &padt, apply_buf);
cannam@127:      else
cannam@127: 	  pln = MKPLAN_HC2C(P, &padt, extra_iter ? apply_extra_iter : apply);
cannam@127: 
cannam@127:      pln->k = ego->k;
cannam@127:      pln->td = 0;
cannam@127:      pln->r = r; pln->rs = X(mkstride)(r, rs);
cannam@127:      pln->m = m; pln->ms = ms;
cannam@127:      pln->v = v; pln->vs = vs;
cannam@127:      pln->slv = ego;
cannam@127:      pln->brs = X(mkstride)(r, 4 * compute_batchsize(r));
cannam@127:      pln->cld0 = cld0;
cannam@127:      pln->cldm = cldm;
cannam@127:      pln->extra_iter = extra_iter;
cannam@127: 
cannam@127:      X(ops_zero)(&pln->super.super.ops);
cannam@127:      X(ops_madd2)(v * (((m - 1) / 2) / e->genus->vl),
cannam@127: 		  &e->ops, &pln->super.super.ops);
cannam@127:      X(ops_madd2)(v, &cld0->ops, &pln->super.super.ops);
cannam@127:      X(ops_madd2)(v, &cldm->ops, &pln->super.super.ops);
cannam@127: 
cannam@127:      if (ego->bufferedp) 
cannam@127: 	  pln->super.super.ops.other += 4 * r * m * v;
cannam@127: 
cannam@127:      return &(pln->super.super);
cannam@127: 
cannam@127:  nada:
cannam@127:      X(plan_destroy_internal)(cld0);
cannam@127:      X(plan_destroy_internal)(cldm);
cannam@127:      return 0;
cannam@127: }
cannam@127: 
cannam@127: static void regone(planner *plnr, khc2c codelet,
cannam@127: 		   const hc2c_desc *desc, 
cannam@127: 		   hc2c_kind hc2ckind, 
cannam@127: 		   int bufferedp)
cannam@127: {
cannam@127:      S *slv = (S *)X(mksolver_hc2c)(sizeof(S), desc->radix, hc2ckind, mkcldw);
cannam@127:      slv->k = codelet;
cannam@127:      slv->desc = desc;
cannam@127:      slv->bufferedp = bufferedp;
cannam@127:      REGISTER_SOLVER(plnr, &(slv->super.super));
cannam@127: }
cannam@127: 
cannam@127: void X(regsolver_hc2c_direct)(planner *plnr, khc2c codelet,
cannam@127: 			      const hc2c_desc *desc,
cannam@127: 			      hc2c_kind hc2ckind)
cannam@127: {
cannam@127:      regone(plnr, codelet, desc, hc2ckind, /* bufferedp */0);
cannam@127:      regone(plnr, codelet, desc, hc2ckind, /* bufferedp */1);
cannam@127: }