annotate fft/fftw/fftw-3.3.4/rdft/rank-geq2-rdft2.c @ 40:223f770b5341 kissfft-double tip

Try a double-precision kissfft
author Chris Cannam
date Wed, 07 Sep 2016 10:40:32 +0100
parents 26056e866c29
children
rev   line source
Chris@19 1 /*
Chris@19 2 * Copyright (c) 2003, 2007-14 Matteo Frigo
Chris@19 3 * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
Chris@19 4 *
Chris@19 5 * This program is free software; you can redistribute it and/or modify
Chris@19 6 * it under the terms of the GNU General Public License as published by
Chris@19 7 * the Free Software Foundation; either version 2 of the License, or
Chris@19 8 * (at your option) any later version.
Chris@19 9 *
Chris@19 10 * This program is distributed in the hope that it will be useful,
Chris@19 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@19 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
Chris@19 13 * GNU General Public License for more details.
Chris@19 14 *
Chris@19 15 * You should have received a copy of the GNU General Public License
Chris@19 16 * along with this program; if not, write to the Free Software
Chris@19 17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
Chris@19 18 *
Chris@19 19 */
Chris@19 20
Chris@19 21
Chris@19 22 /* plans for RDFT2 of rank >= 2 (multidimensional) */
Chris@19 23
Chris@19 24 #include "rdft.h"
Chris@19 25 #include "dft.h"
Chris@19 26
Chris@19 27 typedef struct {
Chris@19 28 solver super;
Chris@19 29 int spltrnk;
Chris@19 30 const int *buddies;
Chris@19 31 int nbuddies;
Chris@19 32 } S;
Chris@19 33
Chris@19 34 typedef struct {
Chris@19 35 plan_dft super;
Chris@19 36 plan *cldr, *cldc;
Chris@19 37 const S *solver;
Chris@19 38 } P;
Chris@19 39
Chris@19 40 static void apply_r2hc(const plan *ego_, R *r0, R *r1, R *cr, R *ci)
Chris@19 41 {
Chris@19 42 const P *ego = (const P *) ego_;
Chris@19 43
Chris@19 44 {
Chris@19 45 plan_rdft2 *cldr = (plan_rdft2 *) ego->cldr;
Chris@19 46 cldr->apply((plan *) cldr, r0, r1, cr, ci);
Chris@19 47 }
Chris@19 48
Chris@19 49 {
Chris@19 50 plan_dft *cldc = (plan_dft *) ego->cldc;
Chris@19 51 cldc->apply((plan *) cldc, cr, ci, cr, ci);
Chris@19 52 }
Chris@19 53 }
Chris@19 54
Chris@19 55 static void apply_hc2r(const plan *ego_, R *r0, R *r1, R *cr, R *ci)
Chris@19 56 {
Chris@19 57 const P *ego = (const P *) ego_;
Chris@19 58
Chris@19 59 {
Chris@19 60 plan_dft *cldc = (plan_dft *) ego->cldc;
Chris@19 61 cldc->apply((plan *) cldc, ci, cr, ci, cr);
Chris@19 62 }
Chris@19 63
Chris@19 64 {
Chris@19 65 plan_rdft2 *cldr = (plan_rdft2 *) ego->cldr;
Chris@19 66 cldr->apply((plan *) cldr, r0, r1, cr, ci);
Chris@19 67 }
Chris@19 68
Chris@19 69 }
Chris@19 70
Chris@19 71 static void awake(plan *ego_, enum wakefulness wakefulness)
Chris@19 72 {
Chris@19 73 P *ego = (P *) ego_;
Chris@19 74 X(plan_awake)(ego->cldr, wakefulness);
Chris@19 75 X(plan_awake)(ego->cldc, wakefulness);
Chris@19 76 }
Chris@19 77
Chris@19 78 static void destroy(plan *ego_)
Chris@19 79 {
Chris@19 80 P *ego = (P *) ego_;
Chris@19 81 X(plan_destroy_internal)(ego->cldr);
Chris@19 82 X(plan_destroy_internal)(ego->cldc);
Chris@19 83 }
Chris@19 84
Chris@19 85 static void print(const plan *ego_, printer *p)
Chris@19 86 {
Chris@19 87 const P *ego = (const P *) ego_;
Chris@19 88 const S *s = ego->solver;
Chris@19 89 p->print(p, "(rdft2-rank>=2/%d%(%p%)%(%p%))",
Chris@19 90 s->spltrnk, ego->cldr, ego->cldc);
Chris@19 91 }
Chris@19 92
Chris@19 93 static int picksplit(const S *ego, const tensor *sz, int *rp)
Chris@19 94 {
Chris@19 95 A(sz->rnk > 1); /* cannot split rnk <= 1 */
Chris@19 96 if (!X(pickdim)(ego->spltrnk, ego->buddies, ego->nbuddies, sz, 1, rp))
Chris@19 97 return 0;
Chris@19 98 *rp += 1; /* convert from dim. index to rank */
Chris@19 99 if (*rp >= sz->rnk) /* split must reduce rank */
Chris@19 100 return 0;
Chris@19 101 return 1;
Chris@19 102 }
Chris@19 103
Chris@19 104 static int applicable0(const solver *ego_, const problem *p_, int *rp,
Chris@19 105 const planner *plnr)
Chris@19 106 {
Chris@19 107 const problem_rdft2 *p = (const problem_rdft2 *) p_;
Chris@19 108 const S *ego = (const S *)ego_;
Chris@19 109 return (1
Chris@19 110 && FINITE_RNK(p->sz->rnk) && FINITE_RNK(p->vecsz->rnk)
Chris@19 111
Chris@19 112 /* FIXME: multidimensional R2HCII ? */
Chris@19 113 && (p->kind == R2HC || p->kind == HC2R)
Chris@19 114
Chris@19 115 && p->sz->rnk >= 2
Chris@19 116 && picksplit(ego, p->sz, rp)
Chris@19 117 && (0
Chris@19 118
Chris@19 119 /* can work out-of-place, but HC2R destroys input */
Chris@19 120 || (p->r0 != p->cr &&
Chris@19 121 (p->kind == R2HC || !NO_DESTROY_INPUTP(plnr)))
Chris@19 122
Chris@19 123 /* FIXME: what are sufficient conditions for inplace? */
Chris@19 124 || (p->r0 == p->cr))
Chris@19 125 );
Chris@19 126 }
Chris@19 127
Chris@19 128 /* TODO: revise this. */
Chris@19 129 static int applicable(const solver *ego_, const problem *p_,
Chris@19 130 const planner *plnr, int *rp)
Chris@19 131 {
Chris@19 132 const S *ego = (const S *)ego_;
Chris@19 133
Chris@19 134 if (!applicable0(ego_, p_, rp, plnr)) return 0;
Chris@19 135
Chris@19 136 if (NO_RANK_SPLITSP(plnr) && (ego->spltrnk != ego->buddies[0]))
Chris@19 137 return 0;
Chris@19 138
Chris@19 139 if (NO_UGLYP(plnr)) {
Chris@19 140 const problem_rdft2 *p = (const problem_rdft2 *) p_;
Chris@19 141
Chris@19 142 /* Heuristic: if the vector stride is greater than the transform
Chris@19 143 size, don't use (prefer to do the vector loop first with a
Chris@19 144 vrank-geq1 plan). */
Chris@19 145 if (p->vecsz->rnk > 0 &&
Chris@19 146 X(tensor_min_stride)(p->vecsz)
Chris@19 147 > X(rdft2_tensor_max_index)(p->sz, p->kind))
Chris@19 148 return 0;
Chris@19 149 }
Chris@19 150
Chris@19 151 return 1;
Chris@19 152 }
Chris@19 153
Chris@19 154 static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
Chris@19 155 {
Chris@19 156 const S *ego = (const S *) ego_;
Chris@19 157 const problem_rdft2 *p;
Chris@19 158 P *pln;
Chris@19 159 plan *cldr = 0, *cldc = 0;
Chris@19 160 tensor *sz1, *sz2, *vecszi, *sz2i;
Chris@19 161 int spltrnk;
Chris@19 162 inplace_kind k;
Chris@19 163 problem *cldp;
Chris@19 164
Chris@19 165 static const plan_adt padt = {
Chris@19 166 X(rdft2_solve), awake, print, destroy
Chris@19 167 };
Chris@19 168
Chris@19 169 if (!applicable(ego_, p_, plnr, &spltrnk))
Chris@19 170 return (plan *) 0;
Chris@19 171
Chris@19 172 p = (const problem_rdft2 *) p_;
Chris@19 173 X(tensor_split)(p->sz, &sz1, spltrnk, &sz2);
Chris@19 174
Chris@19 175 k = p->kind == R2HC ? INPLACE_OS : INPLACE_IS;
Chris@19 176 vecszi = X(tensor_copy_inplace)(p->vecsz, k);
Chris@19 177 sz2i = X(tensor_copy_inplace)(sz2, k);
Chris@19 178
Chris@19 179 /* complex data is ~half of real */
Chris@19 180 sz2i->dims[sz2i->rnk - 1].n = sz2i->dims[sz2i->rnk - 1].n/2 + 1;
Chris@19 181
Chris@19 182 cldr = X(mkplan_d)(plnr,
Chris@19 183 X(mkproblem_rdft2_d)(X(tensor_copy)(sz2),
Chris@19 184 X(tensor_append)(p->vecsz, sz1),
Chris@19 185 p->r0, p->r1,
Chris@19 186 p->cr, p->ci, p->kind));
Chris@19 187 if (!cldr) goto nada;
Chris@19 188
Chris@19 189 if (p->kind == R2HC)
Chris@19 190 cldp = X(mkproblem_dft_d)(X(tensor_copy_inplace)(sz1, k),
Chris@19 191 X(tensor_append)(vecszi, sz2i),
Chris@19 192 p->cr, p->ci, p->cr, p->ci);
Chris@19 193 else /* HC2R must swap re/im parts to get IDFT */
Chris@19 194 cldp = X(mkproblem_dft_d)(X(tensor_copy_inplace)(sz1, k),
Chris@19 195 X(tensor_append)(vecszi, sz2i),
Chris@19 196 p->ci, p->cr, p->ci, p->cr);
Chris@19 197 cldc = X(mkplan_d)(plnr, cldp);
Chris@19 198 if (!cldc) goto nada;
Chris@19 199
Chris@19 200 pln = MKPLAN_RDFT2(P, &padt, p->kind == R2HC ? apply_r2hc : apply_hc2r);
Chris@19 201
Chris@19 202 pln->cldr = cldr;
Chris@19 203 pln->cldc = cldc;
Chris@19 204
Chris@19 205 pln->solver = ego;
Chris@19 206 X(ops_add)(&cldr->ops, &cldc->ops, &pln->super.super.ops);
Chris@19 207
Chris@19 208 X(tensor_destroy4)(sz2i, vecszi, sz2, sz1);
Chris@19 209
Chris@19 210 return &(pln->super.super);
Chris@19 211
Chris@19 212 nada:
Chris@19 213 X(plan_destroy_internal)(cldr);
Chris@19 214 X(plan_destroy_internal)(cldc);
Chris@19 215 X(tensor_destroy4)(sz2i, vecszi, sz2, sz1);
Chris@19 216 return (plan *) 0;
Chris@19 217 }
Chris@19 218
Chris@19 219 static solver *mksolver(int spltrnk, const int *buddies, int nbuddies)
Chris@19 220 {
Chris@19 221 static const solver_adt sadt = { PROBLEM_RDFT2, mkplan, 0 };
Chris@19 222 S *slv = MKSOLVER(S, &sadt);
Chris@19 223 slv->spltrnk = spltrnk;
Chris@19 224 slv->buddies = buddies;
Chris@19 225 slv->nbuddies = nbuddies;
Chris@19 226 return &(slv->super);
Chris@19 227 }
Chris@19 228
Chris@19 229 void X(rdft2_rank_geq2_register)(planner *p)
Chris@19 230 {
Chris@19 231 int i;
Chris@19 232 static const int buddies[] = { 1, 0, -2 };
Chris@19 233
Chris@19 234 const int nbuddies = (int)(sizeof(buddies) / sizeof(buddies[0]));
Chris@19 235
Chris@19 236 for (i = 0; i < nbuddies; ++i)
Chris@19 237 REGISTER_SOLVER(p, mksolver(buddies[i], buddies, nbuddies));
Chris@19 238
Chris@19 239 /* FIXME: Should we try more buddies? See also dft/rank-geq2. */
Chris@19 240 }