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: /* rank-0, vector-rank-3, non-square in-place transposition
Chris@10:    (see rank0.c for square transposition)  */
Chris@10: 
Chris@10: #include "rdft.h"
Chris@10: 
Chris@10: #ifdef HAVE_STRING_H
Chris@10: #include <string.h>		/* for memcpy() */
Chris@10: #endif
Chris@10: 
Chris@10: struct P_s;
Chris@10: 
Chris@10: typedef struct {
Chris@10:      rdftapply apply;
Chris@10:      int (*applicable)(const problem_rdft *p, planner *plnr,
Chris@10: 		       int dim0, int dim1, int dim2, INT *nbuf);
Chris@10:      int (*mkcldrn)(const problem_rdft *p, planner *plnr, struct P_s *ego);
Chris@10:      const char *nam;
Chris@10: } transpose_adt;
Chris@10: 
Chris@10: typedef struct {
Chris@10:      solver super;
Chris@10:      const transpose_adt *adt;
Chris@10: } S;
Chris@10: 
Chris@10: typedef struct P_s {
Chris@10:      plan_rdft super;
Chris@10:      INT n, m, vl; /* transpose n x m matrix of vl-tuples */
Chris@10:      INT nbuf; /* buffer size */
Chris@10:      INT nd, md, d; /* transpose-gcd params */
Chris@10:      INT nc, mc; /* transpose-cut params */
Chris@10:      plan *cld1, *cld2, *cld3; /* children, null if unused */
Chris@10:      const S *slv;
Chris@10: } P;
Chris@10: 
Chris@10: 
Chris@10: /*************************************************************************/
Chris@10: /* some utilities for the solvers */
Chris@10: 
Chris@10: static INT gcd(INT a, INT b)
Chris@10: {
Chris@10:      INT r;
Chris@10:      do {
Chris@10: 	  r = a % b;
Chris@10: 	  a = b;
Chris@10: 	  b = r;
Chris@10:      } while (r != 0);
Chris@10:      
Chris@10:      return a;
Chris@10: }
Chris@10: 
Chris@10: /* whether we can transpose with one of our routines expecting
Chris@10:    contiguous Ntuples */
Chris@10: static int Ntuple_transposable(const iodim *a, const iodim *b, INT vl, INT vs)
Chris@10: {
Chris@10:      return (vs == 1 && b->is == vl && a->os == vl &&
Chris@10: 	     ((a->n == b->n && a->is == b->os
Chris@10: 	       && a->is >= b->n && a->is % vl == 0)
Chris@10: 	      || (a->is == b->n * vl && b->os == a->n * vl)));
Chris@10: }
Chris@10: 
Chris@10: /* check whether a and b correspond to the first and second dimensions
Chris@10:    of a transpose of tuples with vector length = vl, stride = vs. */
Chris@10: static int transposable(const iodim *a, const iodim *b, INT vl, INT vs)
Chris@10: {
Chris@10:      return ((a->n == b->n && a->os == b->is && a->is == b->os)
Chris@10:              || Ntuple_transposable(a, b, vl, vs));
Chris@10: }
Chris@10: 
Chris@10: static int pickdim(const tensor *s, int *pdim0, int *pdim1, int *pdim2)
Chris@10: {
Chris@10:      int dim0, dim1;
Chris@10: 
Chris@10:      for (dim0 = 0; dim0 < s->rnk; ++dim0)
Chris@10:           for (dim1 = 0; dim1 < s->rnk; ++dim1) {
Chris@10: 	       int dim2 = 3 - dim0 - dim1;
Chris@10: 	       if (dim0 == dim1) continue;
Chris@10:                if ((s->rnk == 2 || s->dims[dim2].is == s->dims[dim2].os)
Chris@10: 		   && transposable(s->dims + dim0, s->dims + dim1, 
Chris@10: 				   s->rnk == 2 ? (INT)1 : s->dims[dim2].n,
Chris@10: 				   s->rnk == 2 ? (INT)1 : s->dims[dim2].is)) {
Chris@10:                     *pdim0 = dim0;
Chris@10:                     *pdim1 = dim1;
Chris@10: 		    *pdim2 = dim2;
Chris@10:                     return 1;
Chris@10:                }
Chris@10: 	  }
Chris@10:      return 0;
Chris@10: }
Chris@10: 
Chris@10: #define MINBUFDIV 9 /* min factor by which buffer is smaller than data */
Chris@10: #define MAXBUF 65536 /* maximum non-ugly buffer */
Chris@10: 
Chris@10: /* generic applicability function */
Chris@10: static int applicable(const solver *ego_, const problem *p_, planner *plnr,
Chris@10: 		      int *dim0, int *dim1, int *dim2, INT *nbuf)
Chris@10: {
Chris@10:      const S *ego = (const S *) ego_;
Chris@10:      const problem_rdft *p = (const problem_rdft *) p_;
Chris@10: 
Chris@10:      return (1
Chris@10: 	     && p->I == p->O
Chris@10: 	     && p->sz->rnk == 0
Chris@10: 	     && (p->vecsz->rnk == 2 || p->vecsz->rnk == 3)
Chris@10: 
Chris@10: 	     && pickdim(p->vecsz, dim0, dim1, dim2)
Chris@10: 
Chris@10: 	     /* UGLY if vecloop in wrong order for locality */
Chris@10: 	     && (!NO_UGLYP(plnr) ||
Chris@10: 		 p->vecsz->rnk == 2 ||
Chris@10: 		 X(iabs)(p->vecsz->dims[*dim2].is)
Chris@10: 		 < X(imax)(X(iabs)(p->vecsz->dims[*dim0].is),
Chris@10: 			   X(iabs)(p->vecsz->dims[*dim0].os)))
Chris@10: 
Chris@10: 	     /* SLOW if non-square */
Chris@10: 	     && (!NO_SLOWP(plnr)
Chris@10: 		 || p->vecsz->dims[*dim0].n == p->vecsz->dims[*dim1].n)
Chris@10: 		      
Chris@10: 	     && ego->adt->applicable(p, plnr, *dim0,*dim1,*dim2,nbuf)
Chris@10: 
Chris@10: 	     /* buffers too big are UGLY */
Chris@10: 	     && ((!NO_UGLYP(plnr) && !CONSERVE_MEMORYP(plnr))
Chris@10: 		 || *nbuf <= MAXBUF
Chris@10: 		 || *nbuf * MINBUFDIV <= X(tensor_sz)(p->vecsz))
Chris@10: 	  );
Chris@10: }
Chris@10: 
Chris@10: static void get_transpose_vec(const problem_rdft *p, int dim2, INT *vl,INT *vs)
Chris@10: {
Chris@10:      if (p->vecsz->rnk == 2) {
Chris@10: 	  *vl = 1; *vs = 1;
Chris@10:      }
Chris@10:      else {
Chris@10: 	  *vl = p->vecsz->dims[dim2].n;
Chris@10: 	  *vs = p->vecsz->dims[dim2].is; /* == os */
Chris@10:      }  
Chris@10: }
Chris@10: 
Chris@10: /*************************************************************************/
Chris@10: /* Cache-oblivious in-place transpose of non-square matrices, based 
Chris@10:    on transposes of blocks given by the gcd of the dimensions.
Chris@10: 
Chris@10:    This algorithm is related to algorithm V5 from Murray Dow,
Chris@10:    "Transposing a matrix on a vector computer," Parallel Computing 21
Chris@10:    (12), 1997-2005 (1995), with the modification that we use
Chris@10:    cache-oblivious recursive transpose subroutines (and we derived
Chris@10:    it independently).
Chris@10:    
Chris@10:    For a p x q matrix, this requires scratch space equal to the size
Chris@10:    of the matrix divided by gcd(p,q).  Alternatively, see also the
Chris@10:    "cut" algorithm below, if |p-q| * gcd(p,q) < max(p,q). */
Chris@10: 
Chris@10: static void apply_gcd(const plan *ego_, R *I, R *O)
Chris@10: {
Chris@10:      const P *ego = (const P *) ego_;
Chris@10:      INT n = ego->nd, m = ego->md, d = ego->d;
Chris@10:      INT vl = ego->vl;
Chris@10:      R *buf = (R *)MALLOC(sizeof(R) * ego->nbuf, BUFFERS);
Chris@10:      INT i, num_el = n*m*d*vl;
Chris@10: 
Chris@10:      A(ego->n == n * d && ego->m == m * d);
Chris@10:      UNUSED(O);
Chris@10: 
Chris@10:      /* Transpose the matrix I in-place, where I is an (n*d) x (m*d) matrix
Chris@10: 	of vl-tuples and buf contains n*m*d*vl elements.  
Chris@10: 	
Chris@10: 	In general, to transpose a p x q matrix, you should call this
Chris@10: 	routine with d = gcd(p, q), n = p/d, and m = q/d.  */
Chris@10: 
Chris@10:      A(n > 0 && m > 0 && vl > 0);
Chris@10:      A(d > 1);
Chris@10: 
Chris@10:      /* treat as (d x n) x (d' x m) matrix.  (d' = d) */
Chris@10:      
Chris@10:      /* First, transpose d x (n x d') x m to d x (d' x n) x m,
Chris@10: 	using the buf matrix.  This consists of d transposes
Chris@10: 	of contiguous n x d' matrices of m-tuples. */
Chris@10:      if (n > 1) {
Chris@10: 	  rdftapply cldapply = ((plan_rdft *) ego->cld1)->apply;
Chris@10: 	  for (i = 0; i < d; ++i) {
Chris@10: 	       cldapply(ego->cld1, I + i*num_el, buf);
Chris@10: 	       memcpy(I + i*num_el, buf, num_el*sizeof(R));
Chris@10: 	  }
Chris@10:      }
Chris@10:      
Chris@10:      /* Now, transpose (d x d') x (n x m) to (d' x d) x (n x m), which
Chris@10: 	is a square in-place transpose of n*m-tuples: */
Chris@10:      {
Chris@10: 	  rdftapply cldapply = ((plan_rdft *) ego->cld2)->apply;
Chris@10: 	  cldapply(ego->cld2, I, I);
Chris@10:      }
Chris@10:      
Chris@10:      /* Finally, transpose d' x ((d x n) x m) to d' x (m x (d x n)),
Chris@10: 	using the buf matrix.  This consists of d' transposes
Chris@10: 	of contiguous d*n x m matrices. */
Chris@10:      if (m > 1) {
Chris@10: 	  rdftapply cldapply = ((plan_rdft *) ego->cld3)->apply;
Chris@10: 	  for (i = 0; i < d; ++i) {
Chris@10: 	       cldapply(ego->cld3, I + i*num_el, buf);
Chris@10: 	       memcpy(I + i*num_el, buf, num_el*sizeof(R));
Chris@10: 	  }
Chris@10:      }
Chris@10: 
Chris@10:      X(ifree)(buf);
Chris@10: }
Chris@10: 
Chris@10: static int applicable_gcd(const problem_rdft *p, planner *plnr,
Chris@10: 			  int dim0, int dim1, int dim2, INT *nbuf)
Chris@10: {
Chris@10:      INT n = p->vecsz->dims[dim0].n;
Chris@10:      INT m = p->vecsz->dims[dim1].n;
Chris@10:      INT d, vl, vs;
Chris@10:      get_transpose_vec(p, dim2, &vl, &vs);
Chris@10:      d = gcd(n, m);
Chris@10:      *nbuf = n * (m / d) * vl;
Chris@10:      return (!NO_SLOWP(plnr) /* FIXME: not really SLOW for large 1d ffts */
Chris@10: 	     && n != m
Chris@10: 	     && d > 1
Chris@10: 	     && Ntuple_transposable(p->vecsz->dims + dim0,
Chris@10: 				    p->vecsz->dims + dim1,
Chris@10: 				    vl, vs));
Chris@10: }
Chris@10: 
Chris@10: static int mkcldrn_gcd(const problem_rdft *p, planner *plnr, P *ego)
Chris@10: {
Chris@10:      INT n = ego->nd, m = ego->md, d = ego->d;
Chris@10:      INT vl = ego->vl;
Chris@10:      R *buf = (R *)MALLOC(sizeof(R) * ego->nbuf, BUFFERS);
Chris@10:      INT num_el = n*m*d*vl;
Chris@10: 
Chris@10:      if (n > 1) {
Chris@10: 	  ego->cld1 = X(mkplan_d)(plnr,
Chris@10: 				  X(mkproblem_rdft_0_d)(
Chris@10: 				       X(mktensor_3d)(n, d*m*vl, m*vl,
Chris@10: 						      d, m*vl, n*m*vl,
Chris@10: 						      m*vl, 1, 1),
Chris@10: 				       TAINT(p->I, num_el), buf));
Chris@10: 	  if (!ego->cld1)
Chris@10: 	       goto nada;
Chris@10: 	  X(ops_madd)(d, &ego->cld1->ops, &ego->super.super.ops,
Chris@10: 		      &ego->super.super.ops);
Chris@10: 	  ego->super.super.ops.other += num_el * d * 2;
Chris@10:      }
Chris@10: 
Chris@10:      ego->cld2 = X(mkplan_d)(plnr,
Chris@10: 			     X(mkproblem_rdft_0_d)(
Chris@10: 				  X(mktensor_3d)(d, d*n*m*vl, n*m*vl,
Chris@10: 						 d, n*m*vl, d*n*m*vl,
Chris@10: 						 n*m*vl, 1, 1),
Chris@10: 				  p->I, p->I));
Chris@10:      if (!ego->cld2)
Chris@10: 	  goto nada;
Chris@10:      X(ops_add2)(&ego->cld2->ops, &ego->super.super.ops);
Chris@10: 
Chris@10:      if (m > 1) {
Chris@10: 	  ego->cld3 = X(mkplan_d)(plnr,
Chris@10: 				  X(mkproblem_rdft_0_d)(
Chris@10: 				       X(mktensor_3d)(d*n, m*vl, vl,
Chris@10: 						      m, vl, d*n*vl,
Chris@10: 						      vl, 1, 1),
Chris@10: 				       TAINT(p->I, num_el), buf));
Chris@10: 	  if (!ego->cld3)
Chris@10: 	       goto nada;
Chris@10: 	  X(ops_madd2)(d, &ego->cld3->ops, &ego->super.super.ops);
Chris@10: 	  ego->super.super.ops.other += num_el * d * 2;
Chris@10:      }
Chris@10: 
Chris@10:      X(ifree)(buf);
Chris@10:      return 1;
Chris@10: 
Chris@10:  nada:
Chris@10:      X(ifree)(buf);
Chris@10:      return 0;
Chris@10: }
Chris@10: 
Chris@10: static const transpose_adt adt_gcd =
Chris@10: {
Chris@10:      apply_gcd, applicable_gcd, mkcldrn_gcd,
Chris@10:      "rdft-transpose-gcd"
Chris@10: };
Chris@10: 
Chris@10: /*************************************************************************/
Chris@10: /* Cache-oblivious in-place transpose of non-square n x m matrices,
Chris@10:    based on transposing a sub-matrix first and then transposing the
Chris@10:    remainder(s) with the help of a buffer.  See also transpose-gcd,
Chris@10:    above, if gcd(n,m) is large.
Chris@10: 
Chris@10:    This algorithm is related to algorithm V3 from Murray Dow,
Chris@10:    "Transposing a matrix on a vector computer," Parallel Computing 21
Chris@10:    (12), 1997-2005 (1995), with the modifications that we use
Chris@10:    cache-oblivious recursive transpose subroutines and we have the
Chris@10:    generalization for large |n-m| below.
Chris@10: 
Chris@10:    The best case, and the one described by Dow, is for |n-m| small, in
Chris@10:    which case we transpose a square sub-matrix of size min(n,m),
Chris@10:    handling the remainder via a buffer.  This requires scratch space
Chris@10:    equal to the size of the matrix times |n-m| / max(n,m).
Chris@10: 
Chris@10:    As a generalization when |n-m| is not small, we also support cutting
Chris@10:    *both* dimensions to an nc x mc matrix which is *not* necessarily
Chris@10:    square, but has a large gcd (and can therefore use transpose-gcd).
Chris@10: */
Chris@10: 
Chris@10: static void apply_cut(const plan *ego_, R *I, R *O)
Chris@10: {
Chris@10:      const P *ego = (const P *) ego_;
Chris@10:      INT n = ego->n, m = ego->m, nc = ego->nc, mc = ego->mc, vl = ego->vl;
Chris@10:      INT i;
Chris@10:      R *buf1 = (R *)MALLOC(sizeof(R) * ego->nbuf, BUFFERS);
Chris@10:      UNUSED(O);
Chris@10: 
Chris@10:      if (m > mc) {
Chris@10: 	  ((plan_rdft *) ego->cld1)->apply(ego->cld1, I + mc*vl, buf1);
Chris@10: 	  for (i = 0; i < nc; ++i)
Chris@10: 	       memmove(I + (mc*vl) * i, I + (m*vl) * i, sizeof(R) * (mc*vl));
Chris@10:      }
Chris@10: 
Chris@10:      ((plan_rdft *) ego->cld2)->apply(ego->cld2, I, I); /* nc x mc transpose */
Chris@10:      
Chris@10:      if (n > nc) {
Chris@10: 	  R *buf2 = buf1 + (m-mc)*(nc*vl); /* FIXME: force better alignment? */
Chris@10: 	  memcpy(buf2, I + nc*(m*vl), (n-nc)*(m*vl)*sizeof(R));
Chris@10: 	  for (i = mc-1; i >= 0; --i)
Chris@10: 	       memmove(I + (n*vl) * i, I + (nc*vl) * i, sizeof(R) * (n*vl));
Chris@10: 	  ((plan_rdft *) ego->cld3)->apply(ego->cld3, buf2, I + nc*vl);
Chris@10:      }
Chris@10: 
Chris@10:      if (m > mc) {
Chris@10: 	  if (n > nc)
Chris@10: 	       for (i = mc; i < m; ++i)
Chris@10: 		    memcpy(I + i*(n*vl), buf1 + (i-mc)*(nc*vl),
Chris@10: 			   (nc*vl)*sizeof(R));
Chris@10: 	  else
Chris@10: 	       memcpy(I + mc*(n*vl), buf1, (m-mc)*(n*vl)*sizeof(R));
Chris@10:      }
Chris@10: 
Chris@10:      X(ifree)(buf1);
Chris@10: }
Chris@10: 
Chris@10: /* only cut one dimension if the resulting buffer is small enough */
Chris@10: static int cut1(INT n, INT m, INT vl)
Chris@10: {
Chris@10:      return (X(imax)(n,m) >= X(iabs)(n-m) * MINBUFDIV
Chris@10: 	     || X(imin)(n,m) * X(iabs)(n-m) * vl <= MAXBUF);
Chris@10: }
Chris@10: 
Chris@10: #define CUT_NSRCH 32 /* range of sizes to search for possible cuts */
Chris@10: 
Chris@10: static int applicable_cut(const problem_rdft *p, planner *plnr,
Chris@10: 			  int dim0, int dim1, int dim2, INT *nbuf)
Chris@10: {
Chris@10:      INT n = p->vecsz->dims[dim0].n;
Chris@10:      INT m = p->vecsz->dims[dim1].n;
Chris@10:      INT vl, vs;
Chris@10:      get_transpose_vec(p, dim2, &vl, &vs);
Chris@10:      *nbuf = 0; /* always small enough to be non-UGLY (?) */
Chris@10:      A(MINBUFDIV <= CUT_NSRCH); /* assumed to avoid inf. loops below */
Chris@10:      return (!NO_SLOWP(plnr) /* FIXME: not really SLOW for large 1d ffts? */
Chris@10: 	     && n != m
Chris@10: 	     
Chris@10: 	     /* Don't call transpose-cut recursively (avoid inf. loops):
Chris@10: 	        the non-square sub-transpose produced when !cut1
Chris@10: 	        should always have gcd(n,m) >= min(CUT_NSRCH,n,m),
Chris@10: 	        for which transpose-gcd is applicable */
Chris@10: 	     && (cut1(n, m, vl)
Chris@10: 		 || gcd(n, m) < X(imin)(MINBUFDIV, X(imin)(n,m)))
Chris@10: 
Chris@10: 	     && Ntuple_transposable(p->vecsz->dims + dim0,
Chris@10: 				    p->vecsz->dims + dim1,
Chris@10: 				    vl, vs));
Chris@10: }
Chris@10: 
Chris@10: static int mkcldrn_cut(const problem_rdft *p, planner *plnr, P *ego)
Chris@10: {
Chris@10:      INT n = ego->n, m = ego->m, nc, mc;
Chris@10:      INT vl = ego->vl;
Chris@10:      R *buf;
Chris@10: 
Chris@10:      /* pick the "best" cut */
Chris@10:      if (cut1(n, m, vl)) {
Chris@10: 	  nc = mc = X(imin)(n,m);
Chris@10:      }
Chris@10:      else {
Chris@10: 	  INT dc, ns, ms;
Chris@10: 	  dc = gcd(m, n); nc = n; mc = m;
Chris@10: 	  /* search for cut with largest gcd
Chris@10: 	     (TODO: different optimality criteria? different search range?) */
Chris@10: 	  for (ms = m; ms > 0 && ms > m - CUT_NSRCH; --ms) {
Chris@10: 	       for (ns = n; ns > 0 && ns > n - CUT_NSRCH; --ns) {
Chris@10: 		    INT ds = gcd(ms, ns);
Chris@10: 		    if (ds > dc) {
Chris@10: 			 dc = ds; nc = ns; mc = ms;
Chris@10: 			 if (dc == X(imin)(ns, ms))
Chris@10: 			      break; /* cannot get larger than this */
Chris@10: 		    }
Chris@10: 	       }
Chris@10: 	       if (dc == X(imin)(n, ms))
Chris@10: 		    break; /* cannot get larger than this */
Chris@10: 	  }
Chris@10: 	  A(dc >= X(imin)(CUT_NSRCH, X(imin)(n, m)));
Chris@10:      }
Chris@10:      ego->nc = nc;
Chris@10:      ego->mc = mc;
Chris@10:      ego->nbuf = (m-mc)*(nc*vl) + (n-nc)*(m*vl);
Chris@10: 
Chris@10:      buf = (R *)MALLOC(sizeof(R) * ego->nbuf, BUFFERS);
Chris@10: 
Chris@10:      if (m > mc) {
Chris@10: 	  ego->cld1 = X(mkplan_d)(plnr,
Chris@10: 				  X(mkproblem_rdft_0_d)(
Chris@10: 				       X(mktensor_3d)(nc, m*vl, vl,
Chris@10: 						      m-mc, vl, nc*vl,
Chris@10: 						      vl, 1, 1),
Chris@10: 				       p->I + mc*vl, buf));
Chris@10: 	  if (!ego->cld1)
Chris@10: 	       goto nada;
Chris@10: 	  X(ops_add2)(&ego->cld1->ops, &ego->super.super.ops);
Chris@10:      }
Chris@10: 
Chris@10:      ego->cld2 = X(mkplan_d)(plnr,
Chris@10: 			     X(mkproblem_rdft_0_d)(
Chris@10: 				  X(mktensor_3d)(nc, mc*vl, vl,
Chris@10: 						 mc, vl, nc*vl,
Chris@10: 						 vl, 1, 1),
Chris@10: 				  p->I, p->I));
Chris@10:      if (!ego->cld2)
Chris@10: 	  goto nada;
Chris@10:      X(ops_add2)(&ego->cld2->ops, &ego->super.super.ops);
Chris@10: 
Chris@10:      if (n > nc) {
Chris@10: 	  ego->cld3 = X(mkplan_d)(plnr,
Chris@10: 				  X(mkproblem_rdft_0_d)(
Chris@10: 				       X(mktensor_3d)(n-nc, m*vl, vl,
Chris@10: 						      m, vl, n*vl,
Chris@10: 						      vl, 1, 1),
Chris@10: 				       buf + (m-mc)*(nc*vl), p->I + nc*vl));
Chris@10: 	  if (!ego->cld3)
Chris@10: 	       goto nada;
Chris@10: 	  X(ops_add2)(&ego->cld3->ops, &ego->super.super.ops);
Chris@10:      }
Chris@10: 
Chris@10:      /* memcpy/memmove operations */
Chris@10:      ego->super.super.ops.other += 2 * vl * (nc*mc * ((m > mc) + (n > nc))
Chris@10: 					     + (n-nc)*m + (m-mc)*nc);
Chris@10: 
Chris@10:      X(ifree)(buf);
Chris@10:      return 1;
Chris@10: 
Chris@10:  nada:
Chris@10:      X(ifree)(buf);
Chris@10:      return 0;
Chris@10: }
Chris@10: 
Chris@10: static const transpose_adt adt_cut =
Chris@10: {
Chris@10:      apply_cut, applicable_cut, mkcldrn_cut,
Chris@10:      "rdft-transpose-cut"
Chris@10: };
Chris@10: 
Chris@10: /*************************************************************************/
Chris@10: /* In-place transpose routine from TOMS, which follows the cycles of
Chris@10:    the permutation so that it writes to each location only once.
Chris@10:    Because of cache-line and other issues, however, this routine is
Chris@10:    typically much slower than transpose-gcd or transpose-cut, even
Chris@10:    though the latter do some extra writes.  On the other hand, if the
Chris@10:    vector length is large then the TOMS routine is best.
Chris@10: 
Chris@10:    The TOMS routine also has the advantage of requiring less buffer
Chris@10:    space for the case of gcd(nx,ny) small.  However, in this case it
Chris@10:    has been superseded by the combination of the generalized
Chris@10:    transpose-cut method with the transpose-gcd method, which can
Chris@10:    always transpose with buffers a small fraction of the array size
Chris@10:    regardless of gcd(nx,ny). */
Chris@10: 
Chris@10: /*
Chris@10:  * TOMS Transpose.  Algorithm 513 (Revised version of algorithm 380).
Chris@10:  * 
Chris@10:  * These routines do in-place transposes of arrays.
Chris@10:  * 
Chris@10:  * [ Cate, E.G. and Twigg, D.W., ACM Transactions on Mathematical Software, 
Chris@10:  *   vol. 3, no. 1, 104-110 (1977) ]
Chris@10:  * 
Chris@10:  * C version by Steven G. Johnson (February 1997).
Chris@10:  */
Chris@10: 
Chris@10: /*
Chris@10:  * "a" is a 1D array of length ny*nx*N which constains the nx x ny
Chris@10:  * matrix of N-tuples to be transposed.  "a" is stored in row-major
Chris@10:  * order (last index varies fastest).  move is a 1D array of length
Chris@10:  * move_size used to store information to speed up the process.  The
Chris@10:  * value move_size=(ny+nx)/2 is recommended.  buf should be an array
Chris@10:  * of length 2*N.
Chris@10:  * 
Chris@10:  */
Chris@10: 
Chris@10: static void transpose_toms513(R *a, INT nx, INT ny, INT N,
Chris@10:                               char *move, INT move_size, R *buf)
Chris@10: {
Chris@10:      INT i, im, mn;
Chris@10:      R *b, *c, *d;
Chris@10:      INT ncount;
Chris@10:      INT k;
Chris@10:      
Chris@10:      /* check arguments and initialize: */
Chris@10:      A(ny > 0 && nx > 0 && N > 0 && move_size > 0);
Chris@10:      
Chris@10:      b = buf;
Chris@10:      
Chris@10:      /* Cate & Twigg have a special case for nx == ny, but we don't
Chris@10: 	bother, since we already have special code for this case elsewhere. */
Chris@10: 
Chris@10:      c = buf + N;
Chris@10:      ncount = 2;		/* always at least 2 fixed points */
Chris@10:      k = (mn = ny * nx) - 1;
Chris@10:      
Chris@10:      for (i = 0; i < move_size; ++i)
Chris@10: 	  move[i] = 0;
Chris@10:      
Chris@10:      if (ny >= 3 && nx >= 3)
Chris@10: 	  ncount += gcd(ny - 1, nx - 1) - 1;	/* # fixed points */
Chris@10:      
Chris@10:      i = 1;
Chris@10:      im = ny;
Chris@10:      
Chris@10:      while (1) {
Chris@10: 	  INT i1, i2, i1c, i2c;
Chris@10: 	  INT kmi;
Chris@10: 	  
Chris@10: 	  /** Rearrange the elements of a loop
Chris@10: 	      and its companion loop: **/
Chris@10: 	  
Chris@10: 	  i1 = i;
Chris@10: 	  kmi = k - i;
Chris@10: 	  i1c = kmi;
Chris@10: 	  switch (N) {
Chris@10: 	      case 1:
Chris@10: 		   b[0] = a[i1];
Chris@10: 		   c[0] = a[i1c];
Chris@10: 		   break;
Chris@10: 	      case 2:
Chris@10: 		   b[0] = a[2*i1];
Chris@10: 		   b[1] = a[2*i1+1];
Chris@10: 		   c[0] = a[2*i1c];
Chris@10: 		   c[1] = a[2*i1c+1];
Chris@10: 		   break;
Chris@10: 	      default:
Chris@10: 		   memcpy(b, &a[N * i1], N * sizeof(R));
Chris@10: 		   memcpy(c, &a[N * i1c], N * sizeof(R));
Chris@10: 	  }
Chris@10: 	  while (1) {
Chris@10: 	       i2 = ny * i1 - k * (i1 / nx);
Chris@10: 	       i2c = k - i2;
Chris@10: 	       if (i1 < move_size)
Chris@10: 		    move[i1] = 1;
Chris@10: 	       if (i1c < move_size)
Chris@10: 		    move[i1c] = 1;
Chris@10: 	       ncount += 2;
Chris@10: 	       if (i2 == i)
Chris@10: 		    break;
Chris@10: 	       if (i2 == kmi) {
Chris@10: 		    d = b;
Chris@10: 		    b = c;
Chris@10: 		    c = d;
Chris@10: 		    break;
Chris@10: 	       }
Chris@10: 	       switch (N) {
Chris@10: 		   case 1:
Chris@10: 			a[i1] = a[i2];
Chris@10: 			a[i1c] = a[i2c];
Chris@10: 			break;
Chris@10: 		   case 2:
Chris@10: 			a[2*i1] = a[2*i2];
Chris@10: 			a[2*i1+1] = a[2*i2+1];
Chris@10: 			a[2*i1c] = a[2*i2c];
Chris@10: 			a[2*i1c+1] = a[2*i2c+1];
Chris@10: 			break;
Chris@10: 		   default:
Chris@10: 			memcpy(&a[N * i1], &a[N * i2], 
Chris@10: 			       N * sizeof(R));
Chris@10: 			memcpy(&a[N * i1c], &a[N * i2c], 
Chris@10: 			       N * sizeof(R));
Chris@10: 	       }
Chris@10: 	       i1 = i2;
Chris@10: 	       i1c = i2c;
Chris@10: 	  }
Chris@10: 	  switch (N) {
Chris@10: 	      case 1:
Chris@10: 		   a[i1] = b[0];
Chris@10: 		   a[i1c] = c[0];
Chris@10: 		   break;
Chris@10: 	      case 2:
Chris@10: 		   a[2*i1] = b[0];
Chris@10: 		   a[2*i1+1] = b[1];
Chris@10: 		   a[2*i1c] = c[0];
Chris@10: 		   a[2*i1c+1] = c[1];
Chris@10: 		   break;
Chris@10: 	      default:
Chris@10: 		   memcpy(&a[N * i1], b, N * sizeof(R));
Chris@10: 		   memcpy(&a[N * i1c], c, N * sizeof(R));
Chris@10: 	  }
Chris@10: 	  if (ncount >= mn)
Chris@10: 	       break;	/* we've moved all elements */
Chris@10: 	  
Chris@10: 	  /** Search for loops to rearrange: **/
Chris@10: 	  
Chris@10: 	  while (1) {
Chris@10: 	       INT max = k - i;
Chris@10: 	       ++i;
Chris@10: 	       A(i <= max);
Chris@10: 	       im += ny;
Chris@10: 	       if (im > k)
Chris@10: 		    im -= k;
Chris@10: 	       i2 = im;
Chris@10: 	       if (i == i2)
Chris@10: 		    continue;
Chris@10: 	       if (i >= move_size) {
Chris@10: 		    while (i2 > i && i2 < max) {
Chris@10: 			 i1 = i2;
Chris@10: 			 i2 = ny * i1 - k * (i1 / nx);
Chris@10: 		    }
Chris@10: 		    if (i2 == i)
Chris@10: 			 break;
Chris@10: 	       } else if (!move[i])
Chris@10: 		    break;
Chris@10: 	  }
Chris@10:      }
Chris@10: }
Chris@10: 
Chris@10: static void apply_toms513(const plan *ego_, R *I, R *O)
Chris@10: {
Chris@10:      const P *ego = (const P *) ego_;
Chris@10:      INT n = ego->n, m = ego->m;
Chris@10:      INT vl = ego->vl;
Chris@10:      R *buf = (R *)MALLOC(sizeof(R) * ego->nbuf, BUFFERS);
Chris@10:      UNUSED(O);
Chris@10:      transpose_toms513(I, n, m, vl, (char *) (buf + 2*vl), (n+m)/2, buf);
Chris@10:      X(ifree)(buf);
Chris@10: }
Chris@10: 
Chris@10: static int applicable_toms513(const problem_rdft *p, planner *plnr,
Chris@10: 			   int dim0, int dim1, int dim2, INT *nbuf)
Chris@10: {
Chris@10:      INT n = p->vecsz->dims[dim0].n;
Chris@10:      INT m = p->vecsz->dims[dim1].n;
Chris@10:      INT vl, vs;
Chris@10:      get_transpose_vec(p, dim2, &vl, &vs);
Chris@10:      *nbuf = 2*vl 
Chris@10: 	  + ((n + m) / 2 * sizeof(char) + sizeof(R) - 1) / sizeof(R);
Chris@10:      return (!NO_SLOWP(plnr)
Chris@10: 	     && (vl > 8 || !NO_UGLYP(plnr)) /* UGLY for small vl */
Chris@10: 	     && n != m
Chris@10: 	     && Ntuple_transposable(p->vecsz->dims + dim0,
Chris@10: 				    p->vecsz->dims + dim1,
Chris@10: 				    vl, vs));
Chris@10: }
Chris@10: 
Chris@10: static int mkcldrn_toms513(const problem_rdft *p, planner *plnr, P *ego)
Chris@10: {
Chris@10:      UNUSED(p); UNUSED(plnr);
Chris@10:      /* heuristic so that TOMS algorithm is last resort for small vl */
Chris@10:      ego->super.super.ops.other += ego->n * ego->m * 2 * (ego->vl + 30);
Chris@10:      return 1;
Chris@10: }
Chris@10: 
Chris@10: static const transpose_adt adt_toms513 =
Chris@10: {
Chris@10:      apply_toms513, applicable_toms513, mkcldrn_toms513,
Chris@10:      "rdft-transpose-toms513"
Chris@10: };
Chris@10: 
Chris@10: /*-----------------------------------------------------------------------*/
Chris@10: /*-----------------------------------------------------------------------*/
Chris@10: /* generic stuff: */
Chris@10: 
Chris@10: static void awake(plan *ego_, enum wakefulness wakefulness)
Chris@10: {
Chris@10:      P *ego = (P *) ego_;
Chris@10:      X(plan_awake)(ego->cld1, wakefulness);
Chris@10:      X(plan_awake)(ego->cld2, wakefulness);
Chris@10:      X(plan_awake)(ego->cld3, wakefulness);
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:      p->print(p, "(%s-%Dx%D%v", ego->slv->adt->nam,
Chris@10: 	      ego->n, ego->m, ego->vl);
Chris@10:      if (ego->cld1) p->print(p, "%(%p%)", ego->cld1);
Chris@10:      if (ego->cld2) p->print(p, "%(%p%)", ego->cld2);
Chris@10:      if (ego->cld3) p->print(p, "%(%p%)", ego->cld3);
Chris@10:      p->print(p, ")");
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->cld3);
Chris@10:      X(plan_destroy_internal)(ego->cld2);
Chris@10:      X(plan_destroy_internal)(ego->cld1);
Chris@10: }
Chris@10: 
Chris@10: static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
Chris@10: {
Chris@10:      const S *ego = (const S *) ego_;
Chris@10:      const problem_rdft *p;
Chris@10:      int dim0, dim1, dim2;
Chris@10:      INT nbuf, vs;
Chris@10:      P *pln;
Chris@10: 
Chris@10:      static const plan_adt padt = {
Chris@10: 	  X(rdft_solve), awake, print, destroy
Chris@10:      };
Chris@10: 
Chris@10:      if (!applicable(ego_, p_, plnr, &dim0, &dim1, &dim2, &nbuf))
Chris@10:           return (plan *) 0;
Chris@10: 
Chris@10:      p = (const problem_rdft *) p_;
Chris@10:      pln = MKPLAN_RDFT(P, &padt, ego->adt->apply);
Chris@10: 
Chris@10:      pln->n = p->vecsz->dims[dim0].n;
Chris@10:      pln->m = p->vecsz->dims[dim1].n;
Chris@10:      get_transpose_vec(p, dim2, &pln->vl, &vs);
Chris@10:      pln->nbuf = nbuf;
Chris@10:      pln->d = gcd(pln->n, pln->m);
Chris@10:      pln->nd = pln->n / pln->d;
Chris@10:      pln->md = pln->m / pln->d;
Chris@10:      pln->slv = ego;
Chris@10: 
Chris@10:      X(ops_zero)(&pln->super.super.ops); /* mkcldrn is responsible for ops */
Chris@10: 
Chris@10:      pln->cld1 = pln->cld2 = pln->cld3 = 0;
Chris@10:      if (!ego->adt->mkcldrn(p, plnr, pln)) {
Chris@10: 	  X(plan_destroy_internal)(&(pln->super.super));
Chris@10: 	  return 0;
Chris@10:      }
Chris@10: 
Chris@10:      return &(pln->super.super);
Chris@10: }
Chris@10: 
Chris@10: static solver *mksolver(const transpose_adt *adt)
Chris@10: {
Chris@10:      static const solver_adt sadt = { PROBLEM_RDFT, mkplan, 0 };
Chris@10:      S *slv = MKSOLVER(S, &sadt);
Chris@10:      slv->adt = adt;
Chris@10:      return &(slv->super);
Chris@10: }
Chris@10: 
Chris@10: void X(rdft_vrank3_transpose_register)(planner *p)
Chris@10: {
Chris@10:      unsigned i;
Chris@10:      static const transpose_adt *const adts[] = {
Chris@10: 	  &adt_gcd, &adt_cut,
Chris@10: 	  &adt_toms513
Chris@10:      };
Chris@10:      for (i = 0; i < sizeof(adts) / sizeof(adts[0]); ++i)
Chris@10:           REGISTER_SOLVER(p, mksolver(adts[i]));
Chris@10: }