Mercurial > hg > sv-dependency-builds

annotate src/fftw-3.3.5/rdft/rank0.c @ 42:2cd0e3b3e1fd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Current fftw source
author Chris Cannam
date Tue, 18 Oct 2016 13:40:26 +0100
parents
children
rev   line source
Chris@42 1 /*
Chris@42 2 * Copyright (c) 2003, 2007-14 Matteo Frigo
Chris@42 3 * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
Chris@42 4 *
Chris@42 5 * This program is free software; you can redistribute it and/or modify
Chris@42 6 * it under the terms of the GNU General Public License as published by
Chris@42 7 * the Free Software Foundation; either version 2 of the License, or
Chris@42 8 * (at your option) any later version.
Chris@42 9 *
Chris@42 10 * This program is distributed in the hope that it will be useful,
Chris@42 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@42 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
Chris@42 13 * GNU General Public License for more details.
Chris@42 14 *
Chris@42 15 * You should have received a copy of the GNU General Public License
Chris@42 16 * along with this program; if not, write to the Free Software
Chris@42 17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
Chris@42 18 *
Chris@42 19 */
Chris@42 20
Chris@42 21
Chris@42 22 /* plans for rank-0 RDFTs (copy operations) */
Chris@42 23
Chris@42 24 #include "rdft.h"
Chris@42 25
Chris@42 26 #ifdef HAVE_STRING_H
Chris@42 27 #include <string.h> /* for memcpy() */
Chris@42 28 #endif
Chris@42 29
Chris@42 30 #define MAXRNK 32 /* FIXME: should malloc() */
Chris@42 31
Chris@42 32 typedef struct {
Chris@42 33 plan_rdft super;
Chris@42 34 INT vl;
Chris@42 35 int rnk;
Chris@42 36 iodim d[MAXRNK];
Chris@42 37 const char *nam;
Chris@42 38 } P;
Chris@42 39
Chris@42 40 typedef struct {
Chris@42 41 solver super;
Chris@42 42 rdftapply apply;
Chris@42 43 int (*applicable)(const P *pln, const problem_rdft *p);
Chris@42 44 const char *nam;
Chris@42 45 } S;
Chris@42 46
Chris@42 47 /* copy up to MAXRNK dimensions from problem into plan. If a
Chris@42 48 contiguous dimension exists, save its length in pln->vl */
Chris@42 49 static int fill_iodim(P *pln, const problem_rdft *p)
Chris@42 50 {
Chris@42 51 int i;
Chris@42 52 const tensor *vecsz = p->vecsz;
Chris@42 53
Chris@42 54 pln->vl = 1;
Chris@42 55 pln->rnk = 0;
Chris@42 56 for (i = 0; i < vecsz->rnk; ++i) {
Chris@42 57 /* extract contiguous dimensions */
Chris@42 58 if (pln->vl == 1 &&
Chris@42 59 vecsz->dims[i].is == 1 && vecsz->dims[i].os == 1)
Chris@42 60 pln->vl = vecsz->dims[i].n;
Chris@42 61 else if (pln->rnk == MAXRNK)
Chris@42 62 return 0;
Chris@42 63 else
Chris@42 64 pln->d[pln->rnk++] = vecsz->dims[i];
Chris@42 65 }
Chris@42 66
Chris@42 67 return 1;
Chris@42 68 }
Chris@42 69
Chris@42 70 /* generic higher-rank copy routine, calls cpy2d() to do the real work */
Chris@42 71 static void copy(const iodim *d, int rnk, INT vl,
Chris@42 72 R *I, R *O,
Chris@42 73 cpy2d_func cpy2d)
Chris@42 74 {
Chris@42 75 A(rnk >= 2);
Chris@42 76 if (rnk == 2)
Chris@42 77 cpy2d(I, O, d[0].n, d[0].is, d[0].os, d[1].n, d[1].is, d[1].os, vl);
Chris@42 78 else {
Chris@42 79 INT i;
Chris@42 80 for (i = 0; i < d[0].n; ++i, I += d[0].is, O += d[0].os)
Chris@42 81 copy(d + 1, rnk - 1, vl, I, O, cpy2d);
Chris@42 82 }
Chris@42 83 }
Chris@42 84
Chris@42 85 /* FIXME: should be more general */
Chris@42 86 static int transposep(const P *pln)
Chris@42 87 {
Chris@42 88 int i;
Chris@42 89
Chris@42 90 for (i = 0; i < pln->rnk - 2; ++i)
Chris@42 91 if (pln->d[i].is != pln->d[i].os)
Chris@42 92 return 0;
Chris@42 93
Chris@42 94 return (pln->d[i].n == pln->d[i+1].n &&
Chris@42 95 pln->d[i].is == pln->d[i+1].os &&
Chris@42 96 pln->d[i].os == pln->d[i+1].is);
Chris@42 97 }
Chris@42 98
Chris@42 99 /* generic higher-rank transpose routine, calls transpose2d() to do
Chris@42 100 * the real work */
Chris@42 101 static void transpose(const iodim *d, int rnk, INT vl,
Chris@42 102 R *I,
Chris@42 103 transpose_func transpose2d)
Chris@42 104 {
Chris@42 105 A(rnk >= 2);
Chris@42 106 if (rnk == 2)
Chris@42 107 transpose2d(I, d[0].n, d[0].is, d[0].os, vl);
Chris@42 108 else {
Chris@42 109 INT i;
Chris@42 110 for (i = 0; i < d[0].n; ++i, I += d[0].is)
Chris@42 111 transpose(d + 1, rnk - 1, vl, I, transpose2d);
Chris@42 112 }
Chris@42 113 }
Chris@42 114
Chris@42 115 /**************************************************************/
Chris@42 116 /* rank 0,1,2, out of place, iterative */
Chris@42 117 static void apply_iter(const plan *ego_, R *I, R *O)
Chris@42 118 {
Chris@42 119 const P *ego = (const P *) ego_;
Chris@42 120
Chris@42 121 switch (ego->rnk) {
Chris@42 122 case 0:
Chris@42 123 X(cpy1d)(I, O, ego->vl, 1, 1, 1);
Chris@42 124 break;
Chris@42 125 case 1:
Chris@42 126 X(cpy1d)(I, O,
Chris@42 127 ego->d[0].n, ego->d[0].is, ego->d[0].os,
Chris@42 128 ego->vl);
Chris@42 129 break;
Chris@42 130 default:
Chris@42 131 copy(ego->d, ego->rnk, ego->vl, I, O, X(cpy2d_ci));
Chris@42 132 break;
Chris@42 133 }
Chris@42 134 }
Chris@42 135
Chris@42 136 static int applicable_iter(const P *pln, const problem_rdft *p)
Chris@42 137 {
Chris@42 138 UNUSED(pln);
Chris@42 139 return (p->I != p->O);
Chris@42 140 }
Chris@42 141
Chris@42 142 /**************************************************************/
Chris@42 143 /* out of place, write contiguous output */
Chris@42 144 static void apply_cpy2dco(const plan *ego_, R *I, R *O)
Chris@42 145 {
Chris@42 146 const P *ego = (const P *) ego_;
Chris@42 147 copy(ego->d, ego->rnk, ego->vl, I, O, X(cpy2d_co));
Chris@42 148 }
Chris@42 149
Chris@42 150 static int applicable_cpy2dco(const P *pln, const problem_rdft *p)
Chris@42 151 {
Chris@42 152 int rnk = pln->rnk;
Chris@42 153 return (1
Chris@42 154 && p->I != p->O
Chris@42 155 && rnk >= 2
Chris@42 156
Chris@42 157 /* must not duplicate apply_iter */
Chris@42 158 && (X(iabs)(pln->d[rnk - 2].is) <= X(iabs)(pln->d[rnk - 1].is)
Chris@42 159 ||
Chris@42 160 X(iabs)(pln->d[rnk - 2].os) <= X(iabs)(pln->d[rnk - 1].os))
Chris@42 161 );
Chris@42 162 }
Chris@42 163
Chris@42 164 /**************************************************************/
Chris@42 165 /* out of place, tiled, no buffering */
Chris@42 166 static void apply_tiled(const plan *ego_, R *I, R *O)
Chris@42 167 {
Chris@42 168 const P *ego = (const P *) ego_;
Chris@42 169 copy(ego->d, ego->rnk, ego->vl, I, O, X(cpy2d_tiled));
Chris@42 170 }
Chris@42 171
Chris@42 172 static int applicable_tiled(const P *pln, const problem_rdft *p)
Chris@42 173 {
Chris@42 174 return (1
Chris@42 175 && p->I != p->O
Chris@42 176 && pln->rnk >= 2
Chris@42 177
Chris@42 178 /* somewhat arbitrary */
Chris@42 179 && X(compute_tilesz)(pln->vl, 1) > 4
Chris@42 180 );
Chris@42 181 }
Chris@42 182
Chris@42 183 /**************************************************************/
Chris@42 184 /* out of place, tiled, with buffer */
Chris@42 185 static void apply_tiledbuf(const plan *ego_, R *I, R *O)
Chris@42 186 {
Chris@42 187 const P *ego = (const P *) ego_;
Chris@42 188 copy(ego->d, ego->rnk, ego->vl, I, O, X(cpy2d_tiledbuf));
Chris@42 189 }
Chris@42 190
Chris@42 191 #define applicable_tiledbuf applicable_tiled
Chris@42 192
Chris@42 193 /**************************************************************/
Chris@42 194 /* rank 0, out of place, using memcpy */
Chris@42 195 static void apply_memcpy(const plan *ego_, R *I, R *O)
Chris@42 196 {
Chris@42 197 const P *ego = (const P *) ego_;
Chris@42 198
Chris@42 199 A(ego->rnk == 0);
Chris@42 200 memcpy(O, I, ego->vl * sizeof(R));
Chris@42 201 }
Chris@42 202
Chris@42 203 static int applicable_memcpy(const P *pln, const problem_rdft *p)
Chris@42 204 {
Chris@42 205 return (1
Chris@42 206 && p->I != p->O
Chris@42 207 && pln->rnk == 0
Chris@42 208 && pln->vl > 2 /* do not bother memcpy-ing complex numbers */
Chris@42 209 );
Chris@42 210 }
Chris@42 211
Chris@42 212 /**************************************************************/
Chris@42 213 /* rank > 0 vecloop, out of place, using memcpy (e.g. out-of-place
Chris@42 214 transposes of vl-tuples ... for large vl it should be more
Chris@42 215 efficient to use memcpy than the tiled stuff). */
Chris@42 216
Chris@42 217 static void memcpy_loop(size_t cpysz, int rnk, const iodim *d, R *I, R *O)
Chris@42 218 {
Chris@42 219 INT i, n = d->n, is = d->is, os = d->os;
Chris@42 220 if (rnk == 1)
Chris@42 221 for (i = 0; i < n; ++i, I += is, O += os)
Chris@42 222 memcpy(O, I, cpysz);
Chris@42 223 else {
Chris@42 224 --rnk; ++d;
Chris@42 225 for (i = 0; i < n; ++i, I += is, O += os)
Chris@42 226 memcpy_loop(cpysz, rnk, d, I, O);
Chris@42 227 }
Chris@42 228 }
Chris@42 229
Chris@42 230 static void apply_memcpy_loop(const plan *ego_, R *I, R *O)
Chris@42 231 {
Chris@42 232 const P *ego = (const P *) ego_;
Chris@42 233 memcpy_loop(ego->vl * sizeof(R), ego->rnk, ego->d, I, O);
Chris@42 234 }
Chris@42 235
Chris@42 236 static int applicable_memcpy_loop(const P *pln, const problem_rdft *p)
Chris@42 237 {
Chris@42 238 return (p->I != p->O
Chris@42 239 && pln->rnk > 0
Chris@42 240 && pln->vl > 2 /* do not bother memcpy-ing complex numbers */);
Chris@42 241 }
Chris@42 242
Chris@42 243 /**************************************************************/
Chris@42 244 /* rank 2, in place, square transpose, iterative */
Chris@42 245 static void apply_ip_sq(const plan *ego_, R *I, R *O)
Chris@42 246 {
Chris@42 247 const P *ego = (const P *) ego_;
Chris@42 248 UNUSED(O);
Chris@42 249 transpose(ego->d, ego->rnk, ego->vl, I, X(transpose));
Chris@42 250 }
Chris@42 251
Chris@42 252
Chris@42 253 static int applicable_ip_sq(const P *pln, const problem_rdft *p)
Chris@42 254 {
Chris@42 255 return (1
Chris@42 256 && p->I == p->O
Chris@42 257 && pln->rnk >= 2
Chris@42 258 && transposep(pln));
Chris@42 259 }
Chris@42 260
Chris@42 261 /**************************************************************/
Chris@42 262 /* rank 2, in place, square transpose, tiled */
Chris@42 263 static void apply_ip_sq_tiled(const plan *ego_, R *I, R *O)
Chris@42 264 {
Chris@42 265 const P *ego = (const P *) ego_;
Chris@42 266 UNUSED(O);
Chris@42 267 transpose(ego->d, ego->rnk, ego->vl, I, X(transpose_tiled));
Chris@42 268 }
Chris@42 269
Chris@42 270 static int applicable_ip_sq_tiled(const P *pln, const problem_rdft *p)
Chris@42 271 {
Chris@42 272 return (1
Chris@42 273 && applicable_ip_sq(pln, p)
Chris@42 274
Chris@42 275 /* somewhat arbitrary */
Chris@42 276 && X(compute_tilesz)(pln->vl, 2) > 4
Chris@42 277 );
Chris@42 278 }
Chris@42 279
Chris@42 280 /**************************************************************/
Chris@42 281 /* rank 2, in place, square transpose, tiled, buffered */
Chris@42 282 static void apply_ip_sq_tiledbuf(const plan *ego_, R *I, R *O)
Chris@42 283 {
Chris@42 284 const P *ego = (const P *) ego_;
Chris@42 285 UNUSED(O);
Chris@42 286 transpose(ego->d, ego->rnk, ego->vl, I, X(transpose_tiledbuf));
Chris@42 287 }
Chris@42 288
Chris@42 289 #define applicable_ip_sq_tiledbuf applicable_ip_sq_tiled
Chris@42 290
Chris@42 291 /**************************************************************/
Chris@42 292 static int applicable(const S *ego, const problem *p_)
Chris@42 293 {
Chris@42 294 const problem_rdft *p = (const problem_rdft *) p_;
Chris@42 295 P pln;
Chris@42 296 return (1
Chris@42 297 && p->sz->rnk == 0
Chris@42 298 && FINITE_RNK(p->vecsz->rnk)
Chris@42 299 && fill_iodim(&pln, p)
Chris@42 300 && ego->applicable(&pln, p)
Chris@42 301 );
Chris@42 302 }
Chris@42 303
Chris@42 304 static void print(const plan *ego_, printer *p)
Chris@42 305 {
Chris@42 306 const P *ego = (const P *) ego_;
Chris@42 307 int i;
Chris@42 308 p->print(p, "(%s/%D", ego->nam, ego->vl);
Chris@42 309 for (i = 0; i < ego->rnk; ++i)
Chris@42 310 p->print(p, "%v", ego->d[i].n);
Chris@42 311 p->print(p, ")");
Chris@42 312 }
Chris@42 313
Chris@42 314 static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
Chris@42 315 {
Chris@42 316 const problem_rdft *p;
Chris@42 317 const S *ego = (const S *) ego_;
Chris@42 318 P *pln;
Chris@42 319 int retval;
Chris@42 320
Chris@42 321 static const plan_adt padt = {
Chris@42 322 X(rdft_solve), X(null_awake), print, X(plan_null_destroy)
Chris@42 323 };
Chris@42 324
Chris@42 325 UNUSED(plnr);
Chris@42 326
Chris@42 327 if (!applicable(ego, p_))
Chris@42 328 return (plan *) 0;
Chris@42 329
Chris@42 330 p = (const problem_rdft *) p_;
Chris@42 331 pln = MKPLAN_RDFT(P, &padt, ego->apply);
Chris@42 332
Chris@42 333 retval = fill_iodim(pln, p);
Chris@42 334 (void)retval; /* UNUSED unless DEBUG */
Chris@42 335 A(retval);
Chris@42 336 A(pln->vl > 0); /* because FINITE_RNK(p->vecsz->rnk) holds */
Chris@42 337 pln->nam = ego->nam;
Chris@42 338
Chris@42 339 /* X(tensor_sz)(p->vecsz) loads, X(tensor_sz)(p->vecsz) stores */
Chris@42 340 X(ops_other)(2 * X(tensor_sz)(p->vecsz), &pln->super.super.ops);
Chris@42 341 return &(pln->super.super);
Chris@42 342 }
Chris@42 343
Chris@42 344
Chris@42 345 void X(rdft_rank0_register)(planner *p)
Chris@42 346 {
Chris@42 347 unsigned i;
Chris@42 348 static struct {
Chris@42 349 rdftapply apply;
Chris@42 350 int (*applicable)(const P *, const problem_rdft *);
Chris@42 351 const char *nam;
Chris@42 352 } tab[] = {
Chris@42 353 { apply_memcpy, applicable_memcpy, "rdft-rank0-memcpy" },
Chris@42 354 { apply_memcpy_loop, applicable_memcpy_loop,
Chris@42 355 "rdft-rank0-memcpy-loop" },
Chris@42 356 { apply_iter, applicable_iter, "rdft-rank0-iter-ci" },
Chris@42 357 { apply_cpy2dco, applicable_cpy2dco, "rdft-rank0-iter-co" },
Chris@42 358 { apply_tiled, applicable_tiled, "rdft-rank0-tiled" },
Chris@42 359 { apply_tiledbuf, applicable_tiledbuf, "rdft-rank0-tiledbuf" },
Chris@42 360 { apply_ip_sq, applicable_ip_sq, "rdft-rank0-ip-sq" },
Chris@42 361 {
Chris@42 362 apply_ip_sq_tiled,
Chris@42 363 applicable_ip_sq_tiled,
Chris@42 364 "rdft-rank0-ip-sq-tiled"
Chris@42 365 },
Chris@42 366 {
Chris@42 367 apply_ip_sq_tiledbuf,
Chris@42 368 applicable_ip_sq_tiledbuf,
Chris@42 369 "rdft-rank0-ip-sq-tiledbuf"
Chris@42 370 },
Chris@42 371 };
Chris@42 372
Chris@42 373 for (i = 0; i < sizeof(tab) / sizeof(tab[0]); ++i) {
Chris@42 374 static const solver_adt sadt = { PROBLEM_RDFT, mkplan, 0 };
Chris@42 375 S *slv = MKSOLVER(S, &sadt);
Chris@42 376 slv->apply = tab[i].apply;
Chris@42 377 slv->applicable = tab[i].applicable;
Chris@42 378 slv->nam = tab[i].nam;
Chris@42 379 REGISTER_SOLVER(p, &(slv->super));
Chris@42 380 }
Chris@42 381 }