1 /* test of fftsg3d.c */
2
3 #include <math.h>
4 #include <stdio.h>
5 #include "alloc.h"
6 #define MAX(x,y) ((x) > (y) ? (x) : (y))
7
8 /* random number generator, 0 <= RND < 1 */
9 #define RND(p) ((*(p) = (*(p) * 7141 + 54773) % 259200) * (1.0 / 259200))
10
11
main()12 int main()
13 {
14 void cdft3d(int, int, int, int, double ***, double *, int *, double *);
15 void rdft3d(int, int, int, int, double ***, double *, int *, double *);
16 void ddct3d(int, int, int, int, double ***, double *, int *, double *);
17 void ddst3d(int, int, int, int, double ***, double *, int *, double *);
18 void putdata3d(int n1, int n2, int n3, double ***a);
19 double errorcheck3d(int n1, int n2, int n3, double scale, double ***a);
20 int *ip, n1, n2, n3, n, nt, i, j;
21 double ***a, *w, err;
22
23 printf("data length n1=? (n1 = power of 2) \n");
24 scanf("%d", &n1);
25 printf("data length n2=? (n2 = power of 2) \n");
26 scanf("%d", &n2);
27 printf("data length n3=? (n3 = power of 2) \n");
28 scanf("%d", &n3);
29
30 a = alloc_3d_double(n1, n2, n3);
31 nt = MAX(n1, n2);
32 n = MAX(nt, n3 / 2);
33 ip = alloc_1d_int(2 + (int) sqrt(n + 0.5));
34 n = MAX(nt, n3) * 3 / 2;
35 w = alloc_1d_double(n);
36 ip[0] = 0;
37
38 /* check of CDFT */
39 putdata3d(n1, n2, n3, a);
40 cdft3d(n1, n2, n3, 1, a, NULL, ip, w);
41 cdft3d(n1, n2, n3, -1, a, NULL, ip, w);
42 err = errorcheck3d(n1, n2, n3, 2.0 / n1 / n2 / n3, a);
43 printf("cdft3d err= %g \n", err);
44
45 /* check of RDFT */
46 putdata3d(n1, n2, n3, a);
47 rdft3d(n1, n2, n3, 1, a, NULL, ip, w);
48 rdft3d(n1, n2, n3, -1, a, NULL, ip, w);
49 err = errorcheck3d(n1, n2, n3, 2.0 / n1 / n2 / n3, a);
50 printf("rdft3d err= %g \n", err);
51
52 /* check of DDCT */
53 putdata3d(n1, n2, n3, a);
54 ddct3d(n1, n2, n3, 1, a, NULL, ip, w);
55 ddct3d(n1, n2, n3, -1, a, NULL, ip, w);
56 for (i = 0; i <= n1 - 1; i++) {
57 for (j = 0; j <= n2 - 1; j++) {
58 a[i][j][0] *= 0.5;
59 }
60 for (j = 0; j <= n3 - 1; j++) {
61 a[i][0][j] *= 0.5;
62 }
63 }
64 for (i = 0; i <= n2 - 1; i++) {
65 for (j = 0; j <= n3 - 1; j++) {
66 a[0][i][j] *= 0.5;
67 }
68 }
69 err = errorcheck3d(n1, n2, n3, 8.0 / n1 / n2 / n3, a);
70 printf("ddct3d err= %g \n", err);
71
72 /* check of DDST */
73 putdata3d(n1, n2, n3, a);
74 ddst3d(n1, n2, n3, 1, a, NULL, ip, w);
75 ddst3d(n1, n2, n3, -1, a, NULL, ip, w);
76 for (i = 0; i <= n1 - 1; i++) {
77 for (j = 0; j <= n2 - 1; j++) {
78 a[i][j][0] *= 0.5;
79 }
80 for (j = 0; j <= n3 - 1; j++) {
81 a[i][0][j] *= 0.5;
82 }
83 }
84 for (i = 0; i <= n2 - 1; i++) {
85 for (j = 0; j <= n3 - 1; j++) {
86 a[0][i][j] *= 0.5;
87 }
88 }
89 err = errorcheck3d(n1, n2, n3, 8.0 / n1 / n2 / n3, a);
90 printf("ddst3d err= %g \n", err);
91
92 free_1d_double(w);
93 free_1d_int(ip);
94 free_3d_double(a);
95 return 0;
96 }
97
98
putdata3d(int n1,int n2,int n3,double *** a)99 void putdata3d(int n1, int n2, int n3, double ***a)
100 {
101 int j1, j2, j3, seed = 0;
102
103 for (j1 = 0; j1 <= n1 - 1; j1++) {
104 for (j2 = 0; j2 <= n2 - 1; j2++) {
105 for (j3 = 0; j3 <= n3 - 1; j3++) {
106 a[j1][j2][j3] = RND(&seed);
107 }
108 }
109 }
110 }
111
112
errorcheck3d(int n1,int n2,int n3,double scale,double *** a)113 double errorcheck3d(int n1, int n2, int n3, double scale, double ***a)
114 {
115 int j1, j2, j3, seed = 0;
116 double err = 0, e;
117
118 for (j1 = 0; j1 <= n1 - 1; j1++) {
119 for (j2 = 0; j2 <= n2 - 1; j2++) {
120 for (j3 = 0; j3 <= n3 - 1; j3++) {
121 e = RND(&seed) - a[j1][j2][j3] * scale;
122 err = MAX(err, fabs(e));
123 }
124 }
125 }
126 return err;
127 }
128
129