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