Rev 107 | Details | Compare with Previous | Last modification | View Log | RSS feed
| Rev | Author | Line No. | Line |
|---|---|---|---|
| 2 | pj | 1 | /* |
| 2 | * Copyright (c) 1997-1999 Massachusetts Institute of Technology |
||
| 3 | * |
||
| 4 | * This program is free software; you can redistribute it and/or modify |
||
| 5 | * it under the terms of the GNU General Public License as published by |
||
| 6 | * the Free Software Foundation; either version 2 of the License, or |
||
| 7 | * (at your option) any later version. |
||
| 8 | * |
||
| 9 | * This program is distributed in the hope that it will be useful, |
||
| 10 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
||
| 11 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
||
| 12 | * GNU General Public License for more details. |
||
| 13 | * |
||
| 14 | * You should have received a copy of the GNU General Public License |
||
| 15 | * along with this program; if not, write to the Free Software |
||
| 16 | * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
||
| 17 | * |
||
| 18 | */ |
||
| 19 | |||
| 107 | pj | 20 | /* $Id: rfftwnd.c,v 1.2 2003-03-24 11:14:59 pj Exp $ */ |
| 2 | pj | 21 | |
| 107 | pj | 22 | #include <fftw-int.h> |
| 23 | #include <rfftw.h> |
||
| 2 | pj | 24 | |
| 25 | /********************** prototypes for rexec2 routines **********************/ |
||
| 26 | |||
| 27 | extern void rfftw_real2c_aux(fftw_plan plan, int howmany, |
||
| 28 | fftw_real *in, int istride, int idist, |
||
| 29 | fftw_complex *out, int ostride, int odist, |
||
| 30 | fftw_real *work); |
||
| 31 | extern void rfftw_c2real_aux(fftw_plan plan, int howmany, |
||
| 32 | fftw_complex *in, int istride, int idist, |
||
| 33 | fftw_real *out, int ostride, int odist, |
||
| 34 | fftw_real *work); |
||
| 35 | extern void rfftw_real2c_overlap_aux(fftw_plan plan, int howmany, |
||
| 36 | fftw_real *in, int istride, int idist, |
||
| 37 | fftw_complex *out, int ostride, int odist, |
||
| 38 | fftw_real *work); |
||
| 39 | extern void rfftw_c2real_overlap_aux(fftw_plan plan, int howmany, |
||
| 40 | fftw_complex *in, int istride, int idist, |
||
| 41 | fftw_real *out, int ostride, int odist, |
||
| 42 | fftw_real *work); |
||
| 43 | |||
| 44 | /********************** Initializing the RFFTWND Plan ***********************/ |
||
| 45 | |||
| 46 | /* |
||
| 47 | * Create an fftwnd_plan specialized for specific arrays. (These |
||
| 48 | * arrays are ignored, however, if they are NULL or if the flags |
||
| 49 | * do not include FFTW_MEASURE.) The main advantage of being |
||
| 50 | * provided arrays like this is that we can do runtime timing |
||
| 51 | * measurements of our options, without worrying about allocating |
||
| 52 | * excessive scratch space. |
||
| 53 | */ |
||
| 54 | fftwnd_plan rfftwnd_create_plan_specific(int rank, const int *n, |
||
| 55 | fftw_direction dir, int flags, |
||
| 56 | fftw_real *in, int istride, |
||
| 57 | fftw_real *out, int ostride) |
||
| 58 | { |
||
| 59 | fftwnd_plan p; |
||
| 60 | int i; |
||
| 61 | int rflags = flags & ~FFTW_IN_PLACE; |
||
| 62 | /* note that we always do rfftw transforms out-of-place in rexec2.c */ |
||
| 63 | |||
| 64 | if (flags & FFTW_IN_PLACE) { |
||
| 65 | out = NULL; |
||
| 66 | ostride = istride; |
||
| 67 | } |
||
| 68 | istride = ostride = 1; /* |
||
| 69 | * strides don't work yet, since it is not |
||
| 70 | * clear whether they apply to real |
||
| 71 | * or complex data |
||
| 72 | */ |
||
| 73 | |||
| 74 | if (!(p = fftwnd_create_plan_aux(rank, n, dir, flags))) |
||
| 75 | return 0; |
||
| 76 | |||
| 77 | for (i = 0; i < rank - 1; ++i) |
||
| 78 | p->n_after[i] = (n[rank - 1]/2 + 1) * (p->n_after[i] / n[rank - 1]); |
||
| 79 | if (rank > 0) |
||
| 80 | p->n[rank - 1] = n[rank - 1] / 2 + 1; |
||
| 81 | |||
| 82 | p->plans = fftwnd_new_plan_array(rank); |
||
| 83 | if (rank > 0 && !p->plans) { |
||
| 84 | rfftwnd_destroy_plan(p); |
||
| 85 | return 0; |
||
| 86 | } |
||
| 87 | if (rank > 0) { |
||
| 88 | p->plans[rank - 1] = rfftw_create_plan(n[rank - 1], dir, rflags); |
||
| 89 | if (!p->plans[rank - 1]) { |
||
| 90 | rfftwnd_destroy_plan(p); |
||
| 91 | return 0; |
||
| 92 | } |
||
| 93 | } |
||
| 94 | if (rank > 1) { |
||
| 95 | if (!(flags & FFTW_MEASURE) || in == 0 |
||
| 96 | || (!p->is_in_place && out == 0)) { |
||
| 97 | if (!fftwnd_create_plans_generic(p->plans, rank - 1, n, |
||
| 98 | dir, flags | FFTW_IN_PLACE)) { |
||
| 99 | rfftwnd_destroy_plan(p); |
||
| 100 | return 0; |
||
| 101 | } |
||
| 102 | } else if (dir == FFTW_COMPLEX_TO_REAL || (flags & FFTW_IN_PLACE)) { |
||
| 103 | if (!fftwnd_create_plans_specific(p->plans, rank - 1, n, |
||
| 104 | p->n_after, |
||
| 105 | dir, flags | FFTW_IN_PLACE, |
||
| 106 | (fftw_complex *) in, |
||
| 107 | istride, |
||
| 108 | 0, 0)) { |
||
| 109 | rfftwnd_destroy_plan(p); |
||
| 110 | return 0; |
||
| 111 | } |
||
| 112 | } else { |
||
| 113 | if (!fftwnd_create_plans_specific(p->plans, rank - 1, n, |
||
| 114 | p->n_after, |
||
| 115 | dir, flags | FFTW_IN_PLACE, |
||
| 116 | (fftw_complex *) out, |
||
| 117 | ostride, |
||
| 118 | 0, 0)) { |
||
| 119 | rfftwnd_destroy_plan(p); |
||
| 120 | return 0; |
||
| 121 | } |
||
| 122 | } |
||
| 123 | } |
||
| 124 | p->nbuffers = 0; |
||
| 125 | p->nwork = fftwnd_work_size(rank, p->n, flags | FFTW_IN_PLACE, |
||
| 126 | p->nbuffers + 1); |
||
| 127 | if (p->nwork && !(flags & FFTW_THREADSAFE)) { |
||
| 128 | p->work = (fftw_complex *) fftw_malloc(p->nwork |
||
| 129 | * sizeof(fftw_complex)); |
||
| 130 | if (!p->work) { |
||
| 131 | rfftwnd_destroy_plan(p); |
||
| 132 | return 0; |
||
| 133 | } |
||
| 134 | } |
||
| 135 | return p; |
||
| 136 | } |
||
| 137 | |||
| 138 | fftwnd_plan rfftw2d_create_plan_specific(int nx, int ny, |
||
| 139 | fftw_direction dir, int flags, |
||
| 140 | fftw_real *in, int istride, |
||
| 141 | fftw_real *out, int ostride) |
||
| 142 | { |
||
| 143 | int n[2]; |
||
| 144 | |||
| 145 | n[0] = nx; |
||
| 146 | n[1] = ny; |
||
| 147 | |||
| 148 | return rfftwnd_create_plan_specific(2, n, dir, flags, |
||
| 149 | in, istride, out, ostride); |
||
| 150 | } |
||
| 151 | |||
| 152 | fftwnd_plan rfftw3d_create_plan_specific(int nx, int ny, int nz, |
||
| 153 | fftw_direction dir, int flags, |
||
| 154 | fftw_real *in, int istride, |
||
| 155 | fftw_real *out, int ostride) |
||
| 156 | { |
||
| 157 | int n[3]; |
||
| 158 | |||
| 159 | n[0] = nx; |
||
| 160 | n[1] = ny; |
||
| 161 | n[2] = nz; |
||
| 162 | |||
| 163 | return rfftwnd_create_plan_specific(3, n, dir, flags, |
||
| 164 | in, istride, out, ostride); |
||
| 165 | } |
||
| 166 | |||
| 167 | /* Create a generic fftwnd plan: */ |
||
| 168 | |||
| 169 | fftwnd_plan rfftwnd_create_plan(int rank, const int *n, |
||
| 170 | fftw_direction dir, int flags) |
||
| 171 | { |
||
| 172 | return rfftwnd_create_plan_specific(rank, n, dir, flags, 0, 1, 0, 1); |
||
| 173 | } |
||
| 174 | |||
| 175 | fftwnd_plan rfftw2d_create_plan(int nx, int ny, |
||
| 176 | fftw_direction dir, int flags) |
||
| 177 | { |
||
| 178 | return rfftw2d_create_plan_specific(nx, ny, dir, flags, 0, 1, 0, 1); |
||
| 179 | } |
||
| 180 | |||
| 181 | fftwnd_plan rfftw3d_create_plan(int nx, int ny, int nz, |
||
| 182 | fftw_direction dir, int flags) |
||
| 183 | { |
||
| 184 | return rfftw3d_create_plan_specific(nx, ny, nz, dir, flags, 0, 1, 0, 1); |
||
| 185 | } |
||
| 186 | |||
| 187 | /************************ Freeing the RFFTWND Plan ************************/ |
||
| 188 | |||
| 189 | void rfftwnd_destroy_plan(fftwnd_plan plan) |
||
| 190 | { |
||
| 191 | fftwnd_destroy_plan(plan); |
||
| 192 | } |
||
| 193 | |||
| 194 | /************************ Printing the RFFTWND Plan ************************/ |
||
| 195 | /* |
||
| 196 | void rfftwnd_fprint_plan(FILE *f, fftwnd_plan plan) |
||
| 197 | { |
||
| 198 | fftwnd_fprint_plan(f, plan); |
||
| 199 | } |
||
| 200 | |||
| 201 | void rfftwnd_print_plan(fftwnd_plan plan) |
||
| 202 | { |
||
| 203 | rfftwnd_fprint_plan(stdout, plan); |
||
| 204 | } |
||
| 205 | */ |
||
| 206 | /*********** Computing the N-Dimensional FFT: Auxiliary Routines ************/ |
||
| 207 | |||
| 208 | void rfftwnd_real2c_aux(fftwnd_plan p, int cur_dim, |
||
| 209 | fftw_real *in, int istride, |
||
| 210 | fftw_complex *out, int ostride, |
||
| 211 | fftw_real *work) |
||
| 212 | { |
||
| 213 | int n_after = p->n_after[cur_dim], n = p->n[cur_dim]; |
||
| 214 | |||
| 215 | if (cur_dim == p->rank - 2) { |
||
| 216 | /* just do the last dimension directly: */ |
||
| 217 | if (p->is_in_place) |
||
| 218 | rfftw_real2c_aux(p->plans[p->rank - 1], n, |
||
| 219 | in, istride, (n_after * istride) * 2, |
||
| 220 | out, istride, n_after * istride, |
||
| 221 | work); |
||
| 222 | else |
||
| 223 | rfftw_real2c_aux(p->plans[p->rank - 1], n, |
||
| 224 | in, istride, p->plans[p->rank - 1]->n * istride, |
||
| 225 | out, ostride, n_after * ostride, |
||
| 226 | work); |
||
| 227 | } else { /* we have at least two dimensions to go */ |
||
| 228 | int nr = p->plans[p->rank - 1]->n; |
||
| 229 | int n_after_r = p->is_in_place ? n_after * 2 |
||
| 230 | : nr * (n_after / (nr/2 + 1)); |
||
| 231 | int i; |
||
| 232 | |||
| 233 | /* |
||
| 234 | * process the subsequent dimensions recursively, in hyperslabs, |
||
| 235 | * to get maximum locality: |
||
| 236 | */ |
||
| 237 | for (i = 0; i < n; ++i) |
||
| 238 | rfftwnd_real2c_aux(p, cur_dim + 1, |
||
| 239 | in + i * n_after_r * istride, istride, |
||
| 240 | out + i * n_after * ostride, ostride, work); |
||
| 241 | } |
||
| 242 | |||
| 243 | /* do the current dimension (in-place): */ |
||
| 244 | fftw(p->plans[cur_dim], n_after, |
||
| 245 | out, n_after * ostride, ostride, |
||
| 246 | (fftw_complex *) work, 1, 0); |
||
| 247 | /* I hate this cast */ |
||
| 248 | } |
||
| 249 | |||
| 250 | void rfftwnd_c2real_aux(fftwnd_plan p, int cur_dim, |
||
| 251 | fftw_complex *in, int istride, |
||
| 252 | fftw_real *out, int ostride, |
||
| 253 | fftw_real *work) |
||
| 254 | { |
||
| 255 | int n_after = p->n_after[cur_dim], n = p->n[cur_dim]; |
||
| 256 | |||
| 257 | /* do the current dimension (in-place): */ |
||
| 258 | fftw(p->plans[cur_dim], n_after, |
||
| 259 | in, n_after * istride, istride, |
||
| 260 | (fftw_complex *) work, 1, 0); |
||
| 261 | |||
| 262 | if (cur_dim == p->rank - 2) { |
||
| 263 | /* just do the last dimension directly: */ |
||
| 264 | if (p->is_in_place) |
||
| 265 | rfftw_c2real_aux(p->plans[p->rank - 1], n, |
||
| 266 | in, istride, n_after * istride, |
||
| 267 | out, istride, (n_after * istride) * 2, |
||
| 268 | work); |
||
| 269 | else |
||
| 270 | rfftw_c2real_aux(p->plans[p->rank - 1], n, |
||
| 271 | in, istride, n_after * istride, |
||
| 272 | out, ostride, p->plans[p->rank - 1]->n * ostride, |
||
| 273 | work); |
||
| 274 | } else { /* we have at least two dimensions to go */ |
||
| 275 | int nr = p->plans[p->rank - 1]->n; |
||
| 276 | int n_after_r = p->is_in_place ? n_after * 2 : |
||
| 277 | nr * (n_after / (nr/2 + 1)); |
||
| 278 | int i; |
||
| 279 | |||
| 280 | /* |
||
| 281 | * process the subsequent dimensions recursively, in hyperslabs, |
||
| 282 | * to get maximum locality: |
||
| 283 | */ |
||
| 284 | for (i = 0; i < n; ++i) |
||
| 285 | rfftwnd_c2real_aux(p, cur_dim + 1, |
||
| 286 | in + i * n_after * istride, istride, |
||
| 287 | out + i * n_after_r * ostride, ostride, work); |
||
| 288 | } |
||
| 289 | } |
||
| 290 | |||
| 291 | /* |
||
| 292 | * alternate version of rfftwnd_aux -- this version pushes the howmany |
||
| 293 | * loop down to the leaves of the computation, for greater locality |
||
| 294 | * in cases where dist < stride. It is also required for correctness |
||
| 295 | * if in==out, and we must call a special version of the executor. |
||
| 296 | * Note that work must point to 'howmany' copies of its data |
||
| 297 | * if in == out. |
||
| 298 | */ |
||
| 299 | |||
| 300 | void rfftwnd_real2c_aux_howmany(fftwnd_plan p, int cur_dim, |
||
| 301 | int howmany, |
||
| 302 | fftw_real *in, int istride, int idist, |
||
| 303 | fftw_complex *out, int ostride, int odist, |
||
| 304 | fftw_complex *work) |
||
| 305 | { |
||
| 306 | int n_after = p->n_after[cur_dim], n = p->n[cur_dim]; |
||
| 307 | int k; |
||
| 308 | |||
| 309 | if (cur_dim == p->rank - 2) { |
||
| 310 | /* just do the last dimension directly: */ |
||
| 311 | if (p->is_in_place) |
||
| 312 | for (k = 0; k < n; ++k) |
||
| 313 | rfftw_real2c_overlap_aux(p->plans[p->rank - 1], howmany, |
||
| 314 | in + (k * n_after * istride) * 2, |
||
| 315 | istride, idist, |
||
| 316 | out + (k * n_after * ostride), |
||
| 317 | ostride, odist, |
||
| 318 | (fftw_real *) work); |
||
| 319 | else { |
||
| 320 | int nlast = p->plans[p->rank - 1]->n; |
||
| 321 | for (k = 0; k < n; ++k) |
||
| 322 | rfftw_real2c_aux(p->plans[p->rank - 1], howmany, |
||
| 323 | in + k * nlast * istride, |
||
| 324 | istride, idist, |
||
| 325 | out + k * n_after * ostride, |
||
| 326 | ostride, odist, |
||
| 327 | (fftw_real *) work); |
||
| 328 | } |
||
| 329 | } else { /* we have at least two dimensions to go */ |
||
| 330 | int nr = p->plans[p->rank - 1]->n; |
||
| 331 | int n_after_r = p->is_in_place ? n_after * 2 : |
||
| 332 | nr * (n_after / (nr/2 + 1)); |
||
| 333 | int i; |
||
| 334 | |||
| 335 | /* |
||
| 336 | * process the subsequent dimensions recursively, in hyperslabs, |
||
| 337 | * to get maximum locality: |
||
| 338 | */ |
||
| 339 | for (i = 0; i < n; ++i) |
||
| 340 | rfftwnd_real2c_aux_howmany(p, cur_dim + 1, howmany, |
||
| 341 | in + i * n_after_r * istride, istride, idist, |
||
| 342 | out + i * n_after * ostride, ostride, odist, |
||
| 343 | work); |
||
| 344 | } |
||
| 345 | |||
| 346 | /* do the current dimension (in-place): */ |
||
| 347 | for (k = 0; k < n_after; ++k) |
||
| 348 | fftw(p->plans[cur_dim], howmany, |
||
| 349 | out + k * ostride, n_after * ostride, odist, |
||
| 350 | work, 1, 0); |
||
| 351 | } |
||
| 352 | |||
| 353 | void rfftwnd_c2real_aux_howmany(fftwnd_plan p, int cur_dim, |
||
| 354 | int howmany, |
||
| 355 | fftw_complex *in, int istride, int idist, |
||
| 356 | fftw_real *out, int ostride, int odist, |
||
| 357 | fftw_complex *work) |
||
| 358 | { |
||
| 359 | int n_after = p->n_after[cur_dim], n = p->n[cur_dim]; |
||
| 360 | int k; |
||
| 361 | |||
| 362 | /* do the current dimension (in-place): */ |
||
| 363 | for (k = 0; k < n_after; ++k) |
||
| 364 | fftw(p->plans[cur_dim], howmany, |
||
| 365 | in + k * istride, n_after * istride, idist, |
||
| 366 | work, 1, 0); |
||
| 367 | |||
| 368 | if (cur_dim == p->rank - 2) { |
||
| 369 | /* just do the last dimension directly: */ |
||
| 370 | if (p->is_in_place) |
||
| 371 | for (k = 0; k < n; ++k) |
||
| 372 | rfftw_c2real_overlap_aux(p->plans[p->rank - 1], howmany, |
||
| 373 | in + (k * n_after * istride), |
||
| 374 | istride, idist, |
||
| 375 | out + (k * n_after * ostride) * 2, |
||
| 376 | ostride, odist, |
||
| 377 | (fftw_real *) work); |
||
| 378 | else { |
||
| 379 | int nlast = p->plans[p->rank - 1]->n; |
||
| 380 | for (k = 0; k < n; ++k) |
||
| 381 | rfftw_c2real_aux(p->plans[p->rank - 1], howmany, |
||
| 382 | in + k * n_after * istride, |
||
| 383 | istride, idist, |
||
| 384 | out + k * nlast * ostride, |
||
| 385 | ostride, odist, |
||
| 386 | (fftw_real *) work); |
||
| 387 | } |
||
| 388 | } else { /* we have at least two dimensions to go */ |
||
| 389 | int nr = p->plans[p->rank - 1]->n; |
||
| 390 | int n_after_r = p->is_in_place ? n_after * 2 |
||
| 391 | : nr * (n_after / (nr/2 + 1)); |
||
| 392 | int i; |
||
| 393 | |||
| 394 | /* |
||
| 395 | * process the subsequent dimensions recursively, in hyperslabs, |
||
| 396 | * to get maximum locality: |
||
| 397 | */ |
||
| 398 | for (i = 0; i < n; ++i) |
||
| 399 | rfftwnd_c2real_aux_howmany(p, cur_dim + 1, howmany, |
||
| 400 | in + i * n_after * istride, istride, idist, |
||
| 401 | out + i * n_after_r * ostride, ostride, odist, |
||
| 402 | work); |
||
| 403 | } |
||
| 404 | } |
||
| 405 | |||
| 406 | /********** Computing the N-Dimensional FFT: User-Visible Routines **********/ |
||
| 407 | |||
| 408 | void rfftwnd_real_to_complex(fftwnd_plan p, int howmany, |
||
| 409 | fftw_real *in, int istride, int idist, |
||
| 410 | fftw_complex *out, int ostride, int odist) |
||
| 411 | { |
||
| 412 | fftw_complex *work = p->work; |
||
| 413 | int rank = p->rank; |
||
| 414 | int free_work = 0; |
||
| 415 | |||
| 416 | if (p->dir != FFTW_REAL_TO_COMPLEX) |
||
| 417 | fftw_die("rfftwnd_real_to_complex with complex-to-real plan"); |
||
| 418 | |||
| 419 | #ifdef FFTW_DEBUG |
||
| 420 | if (p->rank > 0 && (p->plans[0]->flags & FFTW_THREADSAFE) |
||
| 421 | && p->nwork && p->work) |
||
| 422 | fftw_die("bug with FFTW_THREADSAFE flag"); |
||
| 423 | #endif |
||
| 424 | |||
| 425 | if (p->is_in_place) { |
||
| 426 | ostride = istride; |
||
| 427 | odist = (idist == 1) ? 1 : (idist / 2); /* ugh */ |
||
| 428 | out = (fftw_complex *) in; |
||
| 429 | if (howmany > 1 && istride > idist && rank > 0) { |
||
| 430 | int new_nwork; |
||
| 431 | |||
| 432 | new_nwork = p->n[rank - 1] * howmany; |
||
| 433 | if (new_nwork > p->nwork) { |
||
| 434 | work = (fftw_complex *) |
||
| 435 | fftw_malloc(sizeof(fftw_complex) * new_nwork); |
||
| 436 | if (!work) |
||
| 437 | fftw_die("error allocating work array"); |
||
| 438 | free_work = 1; |
||
| 439 | } |
||
| 440 | } |
||
| 441 | } |
||
| 442 | if (p->nwork && !work) { |
||
| 443 | work = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * p->nwork); |
||
| 444 | free_work = 1; |
||
| 445 | } |
||
| 446 | switch (rank) { |
||
| 447 | case 0: |
||
| 448 | break; |
||
| 449 | case 1: |
||
| 450 | if (p->is_in_place && howmany > 1 && istride > idist) |
||
| 451 | rfftw_real2c_overlap_aux(p->plans[0], howmany, |
||
| 452 | in, istride, idist, |
||
| 453 | out, ostride, odist, |
||
| 454 | (fftw_real *) work); |
||
| 455 | else |
||
| 456 | rfftw_real2c_aux(p->plans[0], howmany, |
||
| 457 | in, istride, idist, |
||
| 458 | out, ostride, odist, |
||
| 459 | (fftw_real *) work); |
||
| 460 | break; |
||
| 461 | default: /* rank >= 2 */ |
||
| 462 | { |
||
| 463 | if (howmany > 1 && ostride > odist) |
||
| 464 | rfftwnd_real2c_aux_howmany(p, 0, howmany, |
||
| 465 | in, istride, idist, |
||
| 466 | out, ostride, odist, |
||
| 467 | work); |
||
| 468 | else { |
||
| 469 | int i; |
||
| 470 | |||
| 471 | for (i = 0; i < howmany; ++i) |
||
| 472 | rfftwnd_real2c_aux(p, 0, |
||
| 473 | in + i * idist, istride, |
||
| 474 | out + i * odist, ostride, |
||
| 475 | (fftw_real *) work); |
||
| 476 | } |
||
| 477 | } |
||
| 478 | } |
||
| 479 | |||
| 480 | if (free_work) |
||
| 481 | fftw_free(work); |
||
| 482 | } |
||
| 483 | |||
| 484 | void rfftwnd_complex_to_real(fftwnd_plan p, int howmany, |
||
| 485 | fftw_complex *in, int istride, int idist, |
||
| 486 | fftw_real *out, int ostride, int odist) |
||
| 487 | { |
||
| 488 | fftw_complex *work = p->work; |
||
| 489 | int rank = p->rank; |
||
| 490 | int free_work = 0; |
||
| 491 | |||
| 492 | if (p->dir != FFTW_COMPLEX_TO_REAL) |
||
| 493 | fftw_die("rfftwnd_complex_to_real with real-to-complex plan"); |
||
| 494 | |||
| 495 | #ifdef FFTW_DEBUG |
||
| 496 | if (p->rank > 0 && (p->plans[0]->flags & FFTW_THREADSAFE) |
||
| 497 | && p->nwork && p->work) |
||
| 498 | fftw_die("bug with FFTW_THREADSAFE flag"); |
||
| 499 | #endif |
||
| 500 | |||
| 501 | if (p->is_in_place) { |
||
| 502 | ostride = istride; |
||
| 503 | odist = idist; |
||
| 504 | odist = (idist == 1) ? 1 : (idist * 2); /* ugh */ |
||
| 505 | out = (fftw_real *) in; |
||
| 506 | if (howmany > 1 && istride > idist && rank > 0) { |
||
| 507 | int new_nwork = p->n[rank - 1] * howmany; |
||
| 508 | if (new_nwork > p->nwork) { |
||
| 509 | work = (fftw_complex *) |
||
| 510 | fftw_malloc(sizeof(fftw_complex) * new_nwork); |
||
| 511 | if (!work) |
||
| 512 | fftw_die("error allocating work array"); |
||
| 513 | free_work = 1; |
||
| 514 | } |
||
| 515 | } |
||
| 516 | } |
||
| 517 | if (p->nwork && !work) { |
||
| 518 | work = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * p->nwork); |
||
| 519 | free_work = 1; |
||
| 520 | } |
||
| 521 | switch (rank) { |
||
| 522 | case 0: |
||
| 523 | break; |
||
| 524 | case 1: |
||
| 525 | if (p->is_in_place && howmany > 1 && istride > idist) |
||
| 526 | rfftw_c2real_overlap_aux(p->plans[0], howmany, |
||
| 527 | in, istride, idist, |
||
| 528 | out, ostride, odist, |
||
| 529 | (fftw_real *) work); |
||
| 530 | else |
||
| 531 | rfftw_c2real_aux(p->plans[0], howmany, |
||
| 532 | in, istride, idist, |
||
| 533 | out, ostride, odist, |
||
| 534 | (fftw_real *) work); |
||
| 535 | break; |
||
| 536 | default: /* rank >= 2 */ |
||
| 537 | { |
||
| 538 | if (howmany > 1 && ostride > odist) |
||
| 539 | rfftwnd_c2real_aux_howmany(p, 0, howmany, |
||
| 540 | in, istride, idist, |
||
| 541 | out, ostride, odist, |
||
| 542 | work); |
||
| 543 | else { |
||
| 544 | int i; |
||
| 545 | |||
| 546 | for (i = 0; i < howmany; ++i) |
||
| 547 | rfftwnd_c2real_aux(p, 0, |
||
| 548 | in + i * idist, istride, |
||
| 549 | out + i * odist, ostride, |
||
| 550 | (fftw_real *) work); |
||
| 551 | } |
||
| 552 | } |
||
| 553 | } |
||
| 554 | |||
| 555 | if (free_work) |
||
| 556 | fftw_free(work); |
||
| 557 | } |
||
| 558 | |||
| 559 | void rfftwnd_one_real_to_complex(fftwnd_plan p, |
||
| 560 | fftw_real *in, fftw_complex *out) |
||
| 561 | { |
||
| 562 | rfftwnd_real_to_complex(p, 1, in, 1, 1, out, 1, 1); |
||
| 563 | } |
||
| 564 | |||
| 565 | void rfftwnd_one_complex_to_real(fftwnd_plan p, |
||
| 566 | fftw_complex *in, fftw_real *out) |
||
| 567 | { |
||
| 568 | rfftwnd_complex_to_real(p, 1, in, 1, 1, out, 1, 1); |
||
| 569 | } |