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