annotate fft/fftw/fftw-3.3.4/dft/rank-geq2.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 DFT of rank >= 2 (multidimensional) */
Chris@19 23
Chris@19 24 #include "dft.h"
Chris@19 25
Chris@19 26 typedef struct {
Chris@19 27 solver super;
Chris@19 28 int spltrnk;
Chris@19 29 const int *buddies;
Chris@19 30 int nbuddies;
Chris@19 31 } S;
Chris@19 32
Chris@19 33 typedef struct {
Chris@19 34 plan_dft super;
Chris@19 35
Chris@19 36 plan *cld1, *cld2;
Chris@19 37 const S *solver;
Chris@19 38 } P;
Chris@19 39
Chris@19 40 /* Compute multi-dimensional DFT by applying the two cld plans
Chris@19 41 (lower-rnk DFTs). */
Chris@19 42 static void apply(const plan *ego_, R *ri, R *ii, R *ro, R *io)
Chris@19 43 {
Chris@19 44 const P *ego = (const P *) ego_;
Chris@19 45 plan_dft *cld1, *cld2;
Chris@19 46
Chris@19 47 cld1 = (plan_dft *) ego->cld1;
Chris@19 48 cld1->apply(ego->cld1, ri, ii, ro, io);
Chris@19 49
Chris@19 50 cld2 = (plan_dft *) ego->cld2;
Chris@19 51 cld2->apply(ego->cld2, ro, io, ro, io);
Chris@19 52 }
Chris@19 53
Chris@19 54
Chris@19 55 static void awake(plan *ego_, enum wakefulness wakefulness)
Chris@19 56 {
Chris@19 57 P *ego = (P *) ego_;
Chris@19 58 X(plan_awake)(ego->cld1, wakefulness);
Chris@19 59 X(plan_awake)(ego->cld2, wakefulness);
Chris@19 60 }
Chris@19 61
Chris@19 62 static void destroy(plan *ego_)
Chris@19 63 {
Chris@19 64 P *ego = (P *) ego_;
Chris@19 65 X(plan_destroy_internal)(ego->cld2);
Chris@19 66 X(plan_destroy_internal)(ego->cld1);
Chris@19 67 }
Chris@19 68
Chris@19 69 static void print(const plan *ego_, printer *p)
Chris@19 70 {
Chris@19 71 const P *ego = (const P *) ego_;
Chris@19 72 const S *s = ego->solver;
Chris@19 73 p->print(p, "(dft-rank>=2/%d%(%p%)%(%p%))",
Chris@19 74 s->spltrnk, ego->cld1, ego->cld2);
Chris@19 75 }
Chris@19 76
Chris@19 77 static int picksplit(const S *ego, const tensor *sz, int *rp)
Chris@19 78 {
Chris@19 79 A(sz->rnk > 1); /* cannot split rnk <= 1 */
Chris@19 80 if (!X(pickdim)(ego->spltrnk, ego->buddies, ego->nbuddies, sz, 1, rp))
Chris@19 81 return 0;
Chris@19 82 *rp += 1; /* convert from dim. index to rank */
Chris@19 83 if (*rp >= sz->rnk) /* split must reduce rank */
Chris@19 84 return 0;
Chris@19 85 return 1;
Chris@19 86 }
Chris@19 87
Chris@19 88 static int applicable0(const solver *ego_, const problem *p_, int *rp)
Chris@19 89 {
Chris@19 90 const problem_dft *p = (const problem_dft *) p_;
Chris@19 91 const S *ego = (const S *)ego_;
Chris@19 92 return (1
Chris@19 93 && FINITE_RNK(p->sz->rnk) && FINITE_RNK(p->vecsz->rnk)
Chris@19 94 && p->sz->rnk >= 2
Chris@19 95 && picksplit(ego, p->sz, rp)
Chris@19 96 );
Chris@19 97 }
Chris@19 98
Chris@19 99 /* TODO: revise this. */
Chris@19 100 static int applicable(const solver *ego_, const problem *p_,
Chris@19 101 const planner *plnr, int *rp)
Chris@19 102 {
Chris@19 103 const S *ego = (const S *)ego_;
Chris@19 104 const problem_dft *p = (const problem_dft *) p_;
Chris@19 105
Chris@19 106 if (!applicable0(ego_, p_, rp)) return 0;
Chris@19 107
Chris@19 108 if (NO_RANK_SPLITSP(plnr) && (ego->spltrnk != ego->buddies[0])) return 0;
Chris@19 109
Chris@19 110 /* Heuristic: if the vector stride is greater than the transform
Chris@19 111 sz, don't use (prefer to do the vector loop first with a
Chris@19 112 vrank-geq1 plan). */
Chris@19 113 if (NO_UGLYP(plnr))
Chris@19 114 if (p->vecsz->rnk > 0 &&
Chris@19 115 X(tensor_min_stride)(p->vecsz) > X(tensor_max_index)(p->sz))
Chris@19 116 return 0;
Chris@19 117
Chris@19 118 return 1;
Chris@19 119 }
Chris@19 120
Chris@19 121 static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
Chris@19 122 {
Chris@19 123 const S *ego = (const S *) ego_;
Chris@19 124 const problem_dft *p;
Chris@19 125 P *pln;
Chris@19 126 plan *cld1 = 0, *cld2 = 0;
Chris@19 127 tensor *sz1, *sz2, *vecszi, *sz2i;
Chris@19 128 int spltrnk;
Chris@19 129
Chris@19 130 static const plan_adt padt = {
Chris@19 131 X(dft_solve), awake, print, destroy
Chris@19 132 };
Chris@19 133
Chris@19 134 if (!applicable(ego_, p_, plnr, &spltrnk))
Chris@19 135 return (plan *) 0;
Chris@19 136
Chris@19 137 p = (const problem_dft *) p_;
Chris@19 138 X(tensor_split)(p->sz, &sz1, spltrnk, &sz2);
Chris@19 139 vecszi = X(tensor_copy_inplace)(p->vecsz, INPLACE_OS);
Chris@19 140 sz2i = X(tensor_copy_inplace)(sz2, INPLACE_OS);
Chris@19 141
Chris@19 142 cld1 = X(mkplan_d)(plnr,
Chris@19 143 X(mkproblem_dft_d)(X(tensor_copy)(sz2),
Chris@19 144 X(tensor_append)(p->vecsz, sz1),
Chris@19 145 p->ri, p->ii, p->ro, p->io));
Chris@19 146 if (!cld1) goto nada;
Chris@19 147
Chris@19 148 cld2 = X(mkplan_d)(plnr,
Chris@19 149 X(mkproblem_dft_d)(
Chris@19 150 X(tensor_copy_inplace)(sz1, INPLACE_OS),
Chris@19 151 X(tensor_append)(vecszi, sz2i),
Chris@19 152 p->ro, p->io, p->ro, p->io));
Chris@19 153 if (!cld2) goto nada;
Chris@19 154
Chris@19 155 pln = MKPLAN_DFT(P, &padt, apply);
Chris@19 156
Chris@19 157 pln->cld1 = cld1;
Chris@19 158 pln->cld2 = cld2;
Chris@19 159
Chris@19 160 pln->solver = ego;
Chris@19 161 X(ops_add)(&cld1->ops, &cld2->ops, &pln->super.super.ops);
Chris@19 162
Chris@19 163 X(tensor_destroy4)(sz1, sz2, vecszi, sz2i);
Chris@19 164
Chris@19 165 return &(pln->super.super);
Chris@19 166
Chris@19 167 nada:
Chris@19 168 X(plan_destroy_internal)(cld2);
Chris@19 169 X(plan_destroy_internal)(cld1);
Chris@19 170 X(tensor_destroy4)(sz1, sz2, vecszi, sz2i);
Chris@19 171 return (plan *) 0;
Chris@19 172 }
Chris@19 173
Chris@19 174 static solver *mksolver(int spltrnk, const int *buddies, int nbuddies)
Chris@19 175 {
Chris@19 176 static const solver_adt sadt = { PROBLEM_DFT, mkplan, 0 };
Chris@19 177 S *slv = MKSOLVER(S, &sadt);
Chris@19 178 slv->spltrnk = spltrnk;
Chris@19 179 slv->buddies = buddies;
Chris@19 180 slv->nbuddies = nbuddies;
Chris@19 181 return &(slv->super);
Chris@19 182 }
Chris@19 183
Chris@19 184 void X(dft_rank_geq2_register)(planner *p)
Chris@19 185 {
Chris@19 186 int i;
Chris@19 187 static const int buddies[] = { 1, 0, -2 };
Chris@19 188
Chris@19 189 const int nbuddies = (int)(sizeof(buddies) / sizeof(buddies[0]));
Chris@19 190
Chris@19 191 for (i = 0; i < nbuddies; ++i)
Chris@19 192 REGISTER_SOLVER(p, mksolver(buddies[i], buddies, nbuddies));
Chris@19 193
Chris@19 194 /* FIXME:
Chris@19 195
Chris@19 196 Should we try more buddies?
Chris@19 197
Chris@19 198 Another possible variant is to swap cld1 and cld2 (or rather,
Chris@19 199 to swap their problems; they are not interchangeable because
Chris@19 200 cld2 must be in-place). In past versions of FFTW, however, I
Chris@19 201 seem to recall that such rearrangements have made little or no
Chris@19 202 difference.
Chris@19 203 */
Chris@19 204 }