-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathxfft.c
1735 lines (1598 loc) · 52.1 KB
/
xfft.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*
* xfft.c --
*
* Implement eXtended FFT for Yorick, using Swarztrauber FFT routines
* compiled with Yorick, or FFTW (version 2 or 3).
*
*-----------------------------------------------------------------------------
*
* Copyright (C) 2009-2018: Éric Thiébaut <https://github.com/emmt/XFFT>
*
* This program is free software: you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option)
* any later version.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
* more details.
*
* You should have received a copy of the GNU General Public License along with
* this program. If not, see <http://www.gnu.org/licenses/>.
*
*-----------------------------------------------------------------------------
*/
/* RATIONALE AND IMPLEMENTATION NOTES
FFTW assumes row-major order while Yorick is column-major order. Hence the
dimension lists are reversed when creating FFTW plans. For real-to-complex
or complex-to-real transforms, this means that the "half" dimension of
complex Hermitian arrays corresponds to the first Yorick dimension (the
fastest varying one).
For real-to-complex or complex-to-real transforms, in terms of number of
real values (a complex counting as 2 reals), the first dimension of the
complex array is 2*(N1/2 + 1) where N1 is the first dimension of the real
array. Whatever is N1, 2*(N1/2 + 1) is always larger than N1, hence
input and output arrays must be different.
In the benchmarks, out-of-place transforms are slightly faster than
in-place ones. It is expected that, most of the time, the user of the
FFTW-Yorick interface wants to perform out-of-place transforms (preserving
the input array). However, (1) using aligned memory imposes to allocate
internals workspaces and to copy arrays from/to this workspaces; (2) not
all out-of-place transforms (notably the complex-to-real ones) are able to
preserve the input array. Hence, to minimize the number of copies (which
have a time cost for large arrays) and achieve quasi-optimal performances
for most cases. All plans have their private workspace of size suitable to
store any real/complex array with dimensions corresponding to the plan. If
memory alignement is required (option ALIGN=1 in fftw_new), the workspace
is allocated with fftw_malloc and destroyed with fftw_free; otherwise,
p_malloc and p_free are used. If memory alignement is required, in-place
transforms are always used, on entry the input array is copied (and perhaps
converted) into the workspace, then the in-place transform is computed,
then the contents of the workspace is copied into the output array (taking
care of zero padding and data cconversion). If memory alignement is not
required, out-of-place transforms are always used and the internal
workspace may be used to avoid destroying the contents of the input array.
In FFTW, a complex-to-real transform has the following constraints:
- input and output arrays cannot have the same sizes;
- preserving input for out-of-place transforms is not possible;
- out-of-place transform is slightly faster;
Hence the strategy for complex-to-real transforms is:
- FFTW2: use out-of-place transforms;
- FFTW3: if memory alignement is imposed, use in-place transforms (to
limit the number of copies); otherwise, use out-of-place transforms;
- for out-of-place transforms, if input cannot be destroyed, the internal
buffer is used as a scratch array.
*/
#include <string.h>
#include <stdio.h>
#include <yapi.h>
#include <pstdlib.h>
#include <play.h>
#include <ydata.h>
#ifdef YORICK_CHAR_IS_SIGNED
# define BYTE signed char
#else
# define BYTE unsigned char
#endif
/* Define implementation name and perform minimal checks. */
#if defined(CFFT)
# if defined(FFTW3) || defined(FFTW2) || defined(MKL)
# define _XFFT_MULTIPLE
# endif
#elif defined(FFTW3)
# if defined(CFFT) || defined(FFTW2) || defined(MKL)
# define _XFFT_MULTIPLE
# endif
#elif defined(FFTW2)
# if defined(CFFT) || defined(FFTW3) || defined(MKL)
# define _XFFT_MULTIPLE
# endif
#elif defined(MKL)
# if defined(CFFT) || defined(FFTW2) || defined(FFTW3)
# define _XFFT_MULTIPLE
# endif
#endif
#ifdef _XFFT_MULTIPLE
# error multiple implementations
#endif
#ifndef XFFT_impl
# error macro XFFT_impl must be defined
#endif
#ifndef XFFT_IMPL
# error macro XFFT_IMPL must be defined
#endif
#define _XFFT_STRINGIFY(a) #a
#define _XFFT_EXPAND(a) a
#define _XFFT_JOIN2(a1,a2) a1##a2
#define _XFFT_JOIN3(a1,a2,a3) a1##a2##a3
#define _XFFT_JOIN4(a1,a2,a3,a4) a1##a2##a3##a4
#define _XFFT_JOIN5(a1,a2,a3,a4,a5) a1##a2##a3##a4##a5
#define XFFT_STRINGIFY(a) _XFFT_STRINGIFY(a)
#define XFFT_JOIN(a,b) _XFFT_JOIN2(a,b)
#define XFFT_JOIN2(a1,a2) _XFFT_JOIN2(a1,a2)
#define XFFT_JOIN3(a1,a2,a3) _XFFT_JOIN3(a1,a2,a3)
#define XFFT_JOIN4(a1,a2,a3,a4) _XFFT_JOIN4(a1,a2,a3,a4)
#define XFFT_JOIN5(a1,a2,a3,a4,a5) _XFFT_JOIN5(a1,a2,a3,a4,a5)
#define XFFT_IMPL_NAME XFFT_STRINGIFY(XFFT_impl)
#define XFFT_PKG_NAME "xfft_" XFFT_IMPL_NAME
/* Load definitions for FFTW. */
#ifdef FFTW3
# include <fftw3.h>
# define CHOICE(a,b) b
#else
# if defined(FFTW_PREFIX) && (FFTW_PREFIX != 0)
# ifdef USE_THREADS
# include <drfftw_threads.h>
# else
# include <drfftw.h>
# endif
# else
# ifdef USE_THREADS
# include <rfftw_threads.h>
# else
# include <rfftw.h>
# endif
# endif
# ifdef FFTW_ENABLE_FLOAT
# error only double precision real supported
# endif
# define CHOICE(a,b) a
#endif
#define FALSE 0
#define TRUE 1
typedef struct _xform xform_t;
static void setup_dimlist(xform_t *xform, long dims[], int forward);
static void create_plan(xform_t *xform, int forward);
/* Indices of keywords. */
static long align_index = -1L;
static long dims_index = -1L;
static long impl_index = -1L;
static long kind_index = -1L;
static long nevals_index = -1L;
static long nthreads_index = -1L;
static long planning_index = -1L;
static long rdims_index = -1L;
static long real_index = -1L;
static long zdims_index = -1L;
/* Methods. */
static void free_xform(void *);
static void print_xform(void *);
static void eval_xform(void *, int);
static void extract_xform(void *, char *);
static y_userobj_t xform_class = {
/* type_name: */ "xfft_" XFFT_IMPL_NAME,
/* on_free: */ free_xform,
/* on_print: */ print_xform,
/* on_eval: */ eval_xform,
/* on_extract: */ extract_xform,
/* uo_ops: */ (void *)0
};
#define XFFT_ESTIMATE 0
#define XFFT_MEASURE 1
#define XFFT_PATIENT 2
#define XFFT_EXHAUSTIVE 3
/* These constants are defined (without the leading "Y") in "fftw.i". */
#define XFFT_DIRECT 0
#define XFFT_ADJOINT 1
#define XFFT_CONJUGATE_TRANSPOSE XFFT_ADJOINT
#define XFFT_INVERSE 2
#define XFFT_INVERSE_ADJOINT 3
#define XFFT_ADJOINT_INVERSE XFFT_INVERSE_ADJOINT
#define XFFT_INVERSE_CONJUGATE_TRANSPOSE XFFT_INVERSE_ADJOINT
#define XFFT_FORWARD XFFT_DIRECT
#define XFFT_BACKWARD XFFT_ADJOINT
/* The 11 different transforms implemented by FFTW. */
#define XFFT_FFT -1
#define XFFT_R2HC 0
#define XFFT_HC2R 1
#define XFFT_DHT 2
#define XFFT_REDFT00 3
#define XFFT_REDFT01 4
#define XFFT_REDFT10 5
#define XFFT_REDFT11 6
#define XFFT_RODFT00 7
#define XFFT_RODFT01 8
#define XFFT_RODFT10 9
#define XFFT_RODFT11 10
#define XFFT_IS_R2R(kind) (XFFT_R2HC <= (kind) && (kind) <= XFFT_RODFT11)
/* Global kind of transform. */
#define COMPLEX_TO_COMPLEX 0
#define COMPLEX_TO_REAL 1
#define REAL_TO_COMPLEX 2
#define REAL_TO_REAL 4
struct _xform {
double scale; /* scaling factor to normalize the FFT */
long dims[Y_DIMSIZE];
int kinds[Y_DIMSIZE]; /* kinds of real-real transform */
long r_size; /* number of real's in the real array */
long z_size; /* number of real's in the complex array */
long nevals;
void *forward; /* FFTW plan for forward transform */
void *backward; /* FFTW plan for backward transform */
void *ws; /* workspace */
int nthreads;
int planning;
int kind; /* kind of transform, one of: COMPLEX_TO_COMPLEX,
(COMPLEX_TO_REAL|REAL_TO_COMPLEX) or REAL_TO_REAL */
int align;
};
static void free_xform(void *ptr)
{
xform_t *xform = (xform_t *)ptr;
#ifdef FFTW3
if (xform->forward != NULL) {
fftw_destroy_plan(xform->forward);
}
if (xform->backward != NULL) {
fftw_destroy_plan(xform->backward);
}
if (xform->ws != NULL) {
if (xform->align) {
fftw_free(xform->ws);
} else {
p_free(xform->ws);
}
}
#else /* FFTW2 */
if (xform->kind == (COMPLEX_TO_REAL|REAL_TO_COMPLEX)) {
if (xform->forward != NULL) {
rfftwnd_destroy_plan(xform->forward);
}
if (xform->backward != NULL) {
rfftwnd_destroy_plan(xform->backward);
}
} else {
if (xform->forward != NULL) {
fftwnd_destroy_plan(xform->forward);
}
if (xform->backward != NULL) {
fftwnd_destroy_plan(xform->backward);
}
}
if (xform->ws != NULL) {
p_free(xform->ws);
}
#endif /* FFTW2/FFTW3 */
}
static const char*
xform_kind_name(int kind)
{
switch (kind) {
case XFFT_FFT: return "FFT";
case XFFT_R2HC: return "R2HC";
case XFFT_HC2R: return "HC2R";
case XFFT_DHT: return "DHT"; /* discrete Hartley transform */
case XFFT_REDFT00: return "DCT-I";
case XFFT_REDFT01: return "DCT-III";
case XFFT_REDFT10: return "DCT-II";
case XFFT_REDFT11: return "DCT-IV";
case XFFT_RODFT00: return "DST-I";
case XFFT_RODFT01: return "DST-III";
case XFFT_RODFT10: return "DST-II";
case XFFT_RODFT11: return "DST-IV";
default: return "Unknown";
}
}
static long
xform_logical_size(int kind, long n)
{
if (n < 1) {
y_error("invalid size for real-real transform");
return 0L;
}
switch (kind) {
case XFFT_R2HC: return n;
case XFFT_HC2R: return n;
case XFFT_DHT: return n;
case XFFT_REDFT00: return 2*(n - 1);
case XFFT_REDFT01: return 2*n;
case XFFT_REDFT10: return 2*n;
case XFFT_REDFT11: return 2*n;
case XFFT_RODFT00: return 2*(n + 1);
case XFFT_RODFT01: return 2*n;
case XFFT_RODFT10: return 2*n;
case XFFT_RODFT11: return 2*n;
}
y_error("invalid kind of real-real transform");
return 0L;
}
#ifdef FFTW3
static int
xform_forward_kind(int kind)
{
switch (kind) {
case XFFT_R2HC: return FFTW_R2HC;
case XFFT_HC2R: return FFTW_HC2R;
case XFFT_DHT: return FFTW_DHT;
case XFFT_REDFT00: return FFTW_REDFT00;
case XFFT_REDFT01: return FFTW_REDFT01;
case XFFT_REDFT10: return FFTW_REDFT10;
case XFFT_REDFT11: return FFTW_REDFT11;
case XFFT_RODFT00: return FFTW_RODFT00;
case XFFT_RODFT01: return FFTW_RODFT01;
case XFFT_RODFT10: return FFTW_RODFT10;
case XFFT_RODFT11: return FFTW_RODFT11;
}
y_error("invalid kind of real-real transform");
return -1;
}
#endif
#ifdef FFTW3
static int
xform_backward_kind(int kind)
{
switch (kind) {
case XFFT_R2HC: return FFTW_HC2R;
case XFFT_HC2R: return FFTW_R2HC;
case XFFT_DHT: return FFTW_DHT;
case XFFT_REDFT00: return FFTW_REDFT00;
case XFFT_REDFT01: return FFTW_REDFT10;
case XFFT_REDFT10: return FFTW_REDFT01;
case XFFT_REDFT11: return FFTW_REDFT11;
case XFFT_RODFT00: return FFTW_RODFT00;
case XFFT_RODFT01: return FFTW_RODFT10;
case XFFT_RODFT10: return FFTW_RODFT01;
case XFFT_RODFT11: return FFTW_RODFT11;
}
y_error("invalid kind of real-real transform");
return -1;
}
#endif
static void print_xform(void *ptr)
{
xform_t *xform = (xform_t *)ptr;
y_print("xfft operator ("XFFT_IMPL_NAME" implementation) for ", 0);
if (xform->kind != REAL_TO_REAL) {
y_print("FFT", 1);
} else {
int i, kind = xform->kinds[1];
if (xform->kinds[0] > 0) {
for (i = 2; i <= xform->kinds[0]; ++i) {
if (xform->kinds[i] != kind) {
kind = -1;
break;
}
}
}
if (kind != -1) {
y_print(xform_kind_name(kind), 1);
} else {
y_print("[", 0);
for (i = 1; i <= xform->kinds[0]; ++i) {
if (i > 1) {
y_print(",", 0);
}
y_print(xform_kind_name(xform->kinds[i]), 0);
}
y_print("]", 1);
}
}
}
/* Implement the on_extract method to query a member of the object. */
static void extract_xform(void *ptr, char *member)
{
xform_t *xform = (xform_t *)ptr;
long index = yget_global(member, 0);
if (index == dims_index ||
index == rdims_index ||
index == zdims_index) {
long dims_of_dims[2], *xform_dims, *dims, j, ndims;
xform_dims = xform->dims;
ndims = xform_dims[0];
if (ndims < 0) {
ypush_nil();
} else {
dims_of_dims[0] = 1;
dims_of_dims[1] = ndims + 1;
dims = ypush_l(dims_of_dims);
for (j = 0; j <= ndims; ++j) {
dims[j] = xform_dims[j];
}
if (ndims >= 1 && index == zdims_index &&
xform->kind == (COMPLEX_TO_REAL|REAL_TO_COMPLEX)) {
dims[1] = dims[1]/2L + 1L;
}
}
} else if (index == real_index) {
ypush_int(xform->kind == (COMPLEX_TO_REAL|REAL_TO_COMPLEX) ? TRUE : FALSE);
} else if (index == align_index) {
ypush_int(xform->align ? TRUE : FALSE);
} else if (index == planning_index) {
ypush_long(xform->planning);
} else if (index == nthreads_index) {
ypush_long(xform->nthreads);
} else if (index == nevals_index) {
ypush_long(xform->nevals);
} else if (index == impl_index) {
ypush_q(NULL)[0] = p_strcpy(XFFT_IMPL_NAME);
} else if (index == kind_index) {
if (xform->kind != REAL_TO_REAL) {
ypush_int(XFFT_FFT);
} else {
int i, kind = xform->kinds[1];
if (xform->kinds[0] > 0) {
for (i = 2; i <= xform->kinds[0]; ++i) {
if (xform->kinds[i] != kind) {
kind = -1;
break;
}
}
}
if (kind != -1) {
ypush_int(kind);
} else {
int* kinds;
long dims[2];
dims[0] = 1;
dims[1] = xform->kinds[0];
kinds = ypush_i(dims);
for (i = 1; i <= xform->kinds[0]; ++i) {
kinds[i - 1] = xform->kinds[i];
}
}
}
} else {
ypush_nil();
}
}
/*---------------------------------------------------------------------------*/
/* CONVERSION FUNCTIONS */
/* For optimal speed, these function perform, at the same time, copy,
conversion, zero padding and rescaling of the values. */
#if (Y_CHAR != 0)
# error code assumes that Y_CHAR = 0
#endif
#if (Y_SHORT != 1)
# error code assumes that Y_SHORT = 1
#endif
#if (Y_INT != 2)
# error code assumes that Y_INT = 2
#endif
#if (Y_LONG != 3)
# error code assumes that Y_LONG = 3
#endif
#if (Y_FLOAT != 4)
# error code assumes that Y_FLOAT = 4
#endif
#if (Y_DOUBLE != 5)
# error code assumes that Y_DOUBLE = 5
#endif
#if (Y_COMPLEX != 6)
# error code assumes that Y_COMPLEX = 6
#endif
#define IS_INTEGER(id) ((id) >= Y_CHAR && (id) <= Y_LONG)
/* ZPAD_R2H - Zero-pad and convert a real array to form an half-Hermitian
array.
Arguments are:
ID is the data type of the input array;
INP is the address of the input array;
OUT is the address of the output array;
NTOT is the number of elements of the real array (not the half Hermitian one);
N1 is the leading dimension of the real array (not the half Hermitian one);
A is the scale factor (1.0 for no scaling).
The real array is N1-by-N2 (with N2 the product of the trailing dimensions)
while the half-Hermitian array is H1-by-N2 with H1 = (N1/2 + 1). Hence,
the half-Hermitian array has 2*H1*N2 = K1*N2 reals while the real array has
NTOT = N1*N2 reals. If N1 is odd, K1 = 2*H1 = N1+1; otherwise
K1 = 2*H1 = N1+2.
*/
#ifdef FFTW3
#define ZPAD_R2H(id, out, inp, n1, ntot, ptr) \
zpad_r2h[id](out, inp, n1, ntot, ptr)
#define FUNCTION(X, TYPE) \
static void zpad_r2h_##X(double *dst, \
const void *inp, \
const long n1, \
const long ntot, \
const double a) \
{ \
const double zero = 0.0; \
const TYPE *src = (const TYPE *)inp; \
long i, j, k1, n2 = ntot/n1; \
if ((n1 & 1L) != 0L) { \
k1 = n1 + 1; \
if (a != 1.0) { \
for (j = 0; j < n2; ++j, src += n1, dst += k1) { \
for (i = 0; i < n1; ++i) { \
dst[i] = ((double)src[i])*a; \
} \
dst[n1] = zero; \
} \
} else { \
for (j = 0; j < n2; ++j, src += n1, dst += k1) { \
for (i = 0; i < n1; ++i) { \
dst[i] = src[i]; \
} \
dst[n1] = zero; \
} \
} \
} else { \
k1 = n1 + 2; \
if (a != 1.0) { \
for (j = 0; j < n2; ++j, src += n1, dst += k1) { \
for (i = 0; i < n1; ++i) { \
dst[i] = ((double)src[i])*a; \
} \
dst[n1] = zero; \
dst[n1 + 1] = zero; \
} \
} else { \
for (j = 0; j < n2; ++j, src += n1, dst += k1) { \
for (i = 0; i < n1; ++i) { \
dst[i] = src[i]; \
} \
dst[n1] = zero; \
dst[n1 + 1] = zero; \
} \
} \
} \
}
FUNCTION(c, BYTE);
FUNCTION(s, short);
FUNCTION(i, int);
FUNCTION(l, long);
FUNCTION(f, float);
FUNCTION(d, double);
#undef FUNCTION
static void (*zpad_r2h[])(double *, const void *, const long,
const long, const double) = {
zpad_r2h_c, zpad_r2h_s, zpad_r2h_i, zpad_r2h_l, zpad_r2h_f, zpad_r2h_d
};
#endif /* FFTW3 */
/* In-place scaling of an array (assuming SCALE != 1). */
static void scale_d(double *arr, const long ntot, const double scale)
{
long i;
for (i = 0; i < ntot; ++i) {
arr[i] *= scale;
}
}
static void scale_z(double *arr, const long ntot, const double scale)
{
long i;
for (i = 0; i < ntot; ++i) {
arr[2*i] *= scale;
arr[2*i + 1] *= scale;
}
}
/* Out-of-place copy, convert and scale an array of any type to a double
array. */
#define FUNCTION(X,TYPE) \
static void scale_##X##_to_d(double *dst, const void *inp, \
const long ntot, const double scale) \
{ \
const TYPE *src = (const TYPE *)inp; \
long i; \
if (scale != 1.0) { \
for (i = 0; i < ntot; ++i) { \
dst[i] = src[i]*scale; \
} \
} else { \
for (i = 0; i < ntot; ++i) { \
dst[i] = src[i]; \
} \
} \
}
FUNCTION(c, BYTE)
FUNCTION(s, short)
FUNCTION(i, int)
FUNCTION(l, long)
FUNCTION(f, float)
FUNCTION(d, double)
#undef FUNCTION
static void (*scale_x_to_d[])(double *, const void *,
const long, const double) = {
scale_c_to_d, scale_s_to_d, scale_i_to_d, scale_l_to_d,
scale_f_to_d, scale_d_to_d };
#define SCALE_X_TO_D(id, inp, out, ntot, scale) \
scale_x_to_d[id](inp, out, ntot, scale)
/* Out-of-place copy, convert and scale an array of any type to a complex
array. */
#define FUNCTION(X,TYPE) \
static void scale_##X##_to_z(double *dst, const void *inp, \
const long ntot, const double scale) \
{ \
const double zero = 0.0; \
const TYPE *src = (const TYPE *)inp; \
long i; \
if (scale != 1.0) { \
for (i = 0; i < ntot; ++i) { \
dst[2*i] = src[i]*scale; \
dst[2*i + 1] = zero; \
} \
} else { \
for (i = 0; i < ntot; ++i) { \
dst[2*i] = src[i]; \
dst[2*i + 1] = zero; \
} \
} \
}
FUNCTION(c, BYTE)
FUNCTION(s, short)
FUNCTION(i, int)
FUNCTION(l, long)
FUNCTION(f, float)
FUNCTION(d, double)
#undef FUNCTION
static void scale_z_to_z(double *dst, const void *inp,
const long ntot, const double scale)
{
const double *src = (const double *)inp;
long i;
if (scale != 1.0) {
for (i = 0; i < ntot; ++i) {
dst[2*i] = src[2*i]*scale;
dst[2*i + 1] = src[2*i + 1]*scale;
}
} else {
memcpy(dst, src, (2*sizeof(double))*ntot);
}
}
static void (*scale_x_to_z[])(double *, const void *,
const long, const double) = {
scale_c_to_z, scale_s_to_z, scale_i_to_z, scale_l_to_z,
scale_f_to_z, scale_d_to_z, scale_z_to_z };
#define SCALE_X_TO_Z(id, inp, out, ntot, scale) \
scale_x_to_z[id](inp, out, ntot, scale)
/*---------------------------------------------------------------------------*/
/* Macros for calling out-of-place FFT routines. */
#ifdef FFTW3
# define OUT_OF_PLACE_FFT(INP, OUT) \
fftw_execute_dft(plan, INP, OUT)
# define OUT_OF_PLACE_FFT_R2C(INP, OUT) \
fftw_execute_dft_r2c(plan, INP, OUT)
# define OUT_OF_PLACE_FFT_C2R(INP, OUT) \
fftw_execute_dft_c2r(plan, INP, OUT)
#else /* FFTW2 */
# if USE_THREADS
# define OUT_OF_PLACE_CALL(PFX,SFX,INP,OUT) \
if (xform->nthreads > 1) \
PFX##_threads##SFX(xform->nthreads, plan, INP, OUT); \
else \
PFX##SFX(plan, INP, OUT)
# else
# define OUT_OF_PLACE_CALL(PFX,SFX,INP,OUT) \
PFX##SFX(plan, INP, OUT)
# endif
# define OUT_OF_PLACE_FFT(INP, OUT) \
OUT_OF_PLACE_CALL(fftwnd,_one,INP,OUT)
# define OUT_OF_PLACE_FFT_R2C(INP, OUT) \
OUT_OF_PLACE_CALL(rfftwnd,_one_real_to_complex,INP,OUT)
# define OUT_OF_PLACE_FFT_C2R(INP, OUT) \
OUT_OF_PLACE_CALL(rfftwnd,_one_complex_to_real,INP,OUT)
#endif /* FFTW3 or FFTW2 */
/* iarg = argc - k is stack index of k-th argument of the builtin function */
/*
* rank = 0, scalar array, FFT is a no-op
* rank = 1
*/
static void eval_xform(void *ptr, int argc)
{
double scale;
xform_t *xform = (xform_t *)ptr;
void *inp, *out, *tmp, *plan;
long ntot, index;
#ifdef FFTW3
long inp_dim1, out_dim1;
int rank;
#endif
int arg_type, xform_kind, forward, overwrite, rescale, job, scratch;
long dims[Y_DIMSIZE];
/* Get the direction of the transform and check number of arguments. */
if (argc == 2) {
arg_type = yarg_typeid(0);
if (IS_INTEGER(arg_type) && yarg_rank(0) == 0) {
job = ygets_i(0);
} else if (arg_type == Y_VOID) {
job = 0;
} else {
job = -1;
}
switch (job) {
case XFFT_DIRECT:
forward = TRUE;
rescale = FALSE;
break;
case XFFT_ADJOINT:
forward = FALSE;
rescale = FALSE;
break;
case XFFT_INVERSE:
forward = FALSE;
rescale = TRUE;
break;
case XFFT_INVERSE_ADJOINT:
forward = TRUE;
rescale = TRUE;
break;
default:
y_error("bad job");
return; /* to avoid compiler warnings */
}
yarg_drop(1);
} else {
/* Default is same as DIRECT. */
if (argc != 1) {
y_error("syntax: op(a) or op(a, job) with op the XFFT operator");
}
forward = TRUE;
rescale = FALSE;
}
/* Figure out whether or not to perform in-place operation and, if needed,
set INDEX to save result in a variable. FIXME: The rationale of
yarg_scratch is strange: it returns 1 (hence true) when the argument is a
variable reference; this is why, I combine the result of yarg_scratch and
yget_ref to properly figure out whether an argument can be used
in-place. */
index = yget_ref(0);
scratch = (index < 0 ? yarg_scratch(0) : FALSE);
if (yarg_subroutine()) {
/* FIXME: it should be possible to perform in-place operation for other
kind of arguments such as hash table members. */
if (index < 0) {
y_error("when called as a subroutine, argument must be "
"a simple variable, not a temporary expression");
}
overwrite = TRUE;
} else {
overwrite = (scratch != 0);
index = -1L; /* avoids setting global symbol at the end */
}
/* Check data type of input array. */
inp = ygeta_any(0, &ntot, dims, &arg_type);
if (arg_type < 0 || arg_type > Y_COMPLEX) {
y_error("invalid non-numerical data type");
}
if (xform->kind == (REAL_TO_COMPLEX|COMPLEX_TO_REAL)) {
if (forward) {
if (arg_type == Y_COMPLEX) {
y_error("invalid complex input for forward real-complex transform");
}
xform_kind = REAL_TO_COMPLEX;
} else {
xform_kind = COMPLEX_TO_REAL;
}
} else {
xform_kind = xform->kind;
if (xform_kind == REAL_TO_REAL && arg_type == Y_COMPLEX) {
y_error("invalid complex input for real-real transform");
}
}
/* Check dimension list of the input array and get dimension list of the
result. Must be done prior to workspace allocation and plan creation
because this may set the dimension list of the transform. */
#ifdef FFTW3
rank = dims[0];
inp_dim1 = (rank > 0 ? dims[1] : 0); /* save 1st dimension of input */
#endif
setup_dimlist(xform, dims, forward);
#ifdef FFTW3
out_dim1 = (rank > 0 ? dims[1] : 0); /* save 1st dimension of output */
#endif
if (ntot == 1) {
/* Transform of a scalar is the identity so, at most, it is just a matter
of converting the data type of the input. */
if (xform_kind == COMPLEX_TO_REAL || xform_kind == REAL_TO_REAL) {
out = ygeta_d(0, NULL, NULL);
} else {
out = ygeta_z(0, NULL, NULL);
}
goto done;
}
/* Choose the plan, allocate it if necessary. */
if (forward) {
if (xform->forward == NULL) create_plan(xform, TRUE);
plan = xform->forward;
} else {
if (xform->backward == NULL) create_plan(xform, FALSE);
plan = xform->backward;
}
scale = (rescale ? xform->scale : 1.0);
if (xform_kind == COMPLEX_TO_COMPLEX) {
/**************************/
/* COMPLEX-TO-COMPLEX FFT */
/**************************/
#ifdef FFTW3
if (xform->align) {
/* Perform in-place transform with internal workspace. */
tmp = xform->ws;
SCALE_X_TO_Z(arg_type, tmp, inp, ntot, scale);
fftw_execute(plan);
if (overwrite && arg_type == Y_COMPLEX) {
out = inp;
index = -1L; /* no needs to set global symbol at the end */
} else {
yarg_drop(1);
out = ypush_z(dims);
}
memcpy(out, tmp, (2*sizeof(double))*ntot);
goto done;
}
#endif
/* Perform out-of-place transform with destroyable input. */
if (overwrite && arg_type == Y_COMPLEX) {
tmp = inp;
if (scale != 1.0) {
scale_z(tmp, ntot, scale);
}
} else {
tmp = xform->ws;
SCALE_X_TO_Z(arg_type, tmp, inp, ntot, scale);
yarg_drop(1);
}
out = ypush_z(dims);
OUT_OF_PLACE_FFT(tmp, out);
} else if (xform_kind == REAL_TO_COMPLEX) {
/***********************/
/* REAL-TO-COMPLEX FFT */
/***********************/
#ifdef FFTW3
if (xform->align) {
tmp = xform->ws;
ZPAD_R2H(arg_type, tmp, inp, inp_dim1, ntot, scale);
yarg_drop(1);
fftw_execute(plan); /* FIXME: fftw_execute_dft_r2c(plan, tmp, tmp); */
out = ypush_z(dims);
memcpy(out, tmp, xform->z_size*sizeof(double));
goto done;
}
#endif
if (overwrite && arg_type == Y_DOUBLE) {
tmp = inp;
if (scale != 1.0) {
scale_d(tmp, ntot, scale);
}
} else {
tmp = xform->ws;
SCALE_X_TO_D(arg_type, tmp, inp, ntot, scale);
yarg_drop(1);
}
out = ypush_z(dims);
OUT_OF_PLACE_FFT_R2C(tmp, out);
} else if (xform_kind == COMPLEX_TO_REAL) {
/***********************/
/* COMPLEX-TO-REAL FFT */
/***********************/
#ifdef FFTW3
if (xform->align) {
/* Perform in-place transform, then perform a packed copy of the
half-Hermitian result into a real array. */
double *dst;
const double *src;
long i, j, k1, n1, n2;
tmp = xform->ws;
SCALE_X_TO_Z(arg_type, tmp, inp, ntot, scale);
fftw_execute(plan); /* FIXME: fftw_execute_dft_c2r(plan, tmp, tmp); */
src = (const double *)tmp;
dst = (double *)ypush_d(dims);
n1 = out_dim1;
n2 = ntot/inp_dim1;
k1 = 2L*inp_dim1;
for (j = 0L; j < n2; ++j, src += k1, dst += n1) {
for (i = 0L; i < n1; ++i) {
dst[i] = src[i];
}
}
goto done;
}
#endif
/* Perform out-of-place transform. */
if (overwrite && arg_type == Y_COMPLEX) {
tmp = inp;
if (scale != 1.0) {
scale_z(tmp, ntot, scale);
}
} else {
tmp = xform->ws;
SCALE_X_TO_Z(arg_type, tmp, inp, ntot, scale);
yarg_drop(1);
}
out = ypush_d(dims);
OUT_OF_PLACE_FFT_C2R(tmp, out);
} else if (xform_kind == REAL_TO_REAL) {
/**************************/
/* REAL-TO-REAL TRANSFORM */
/**************************/
#ifdef FFTW3
if (xform->align) {
/* Perform in-place transform. */
SCALE_X_TO_D(arg_type, xform->ws, inp, ntot, scale);
fftw_execute(plan);
(void)memcpy(ypush_d(dims), xform->ws, ntot*sizeof(double));
} else {
/* Perform out-of-place transform. */
if (overwrite && arg_type == Y_DOUBLE) {
tmp = inp;
if (scale != 1.0) {
scale_d(tmp, ntot, scale);
}