NFFT  3.3.0
fastsum_benchomp_createdataset.c
1 /*
2  * Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
3  *
4  * This program is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU General Public License as published by the Free Software
6  * Foundation; either version 2 of the License, or (at your option) any later
7  * version.
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU General Public License along with
15  * this program; if not, write to the Free Software Foundation, Inc., 51
16  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18 
19 /* $Id: simple_test.c 3372 2009-10-21 06:04:05Z skunis $ */
20 
21 #include <stdio.h>
22 #include <math.h>
23 #include <string.h>
24 #include <stdlib.h>
25 #include <complex.h>
26 
27 #include "config.h"
28 
29 #include "nfft3.h"
30 #include "infft.h"
31 
32 void fastsum_benchomp_createdataset(unsigned int d, int L, int M)
33 {
34  int t, j, k;
35  R *x;
36  R *y;
37  C *alpha;
38 
39  x = (R*) NFFT(malloc)((size_t)(d * L) * sizeof(R));
40  y = (R*) NFFT(malloc)((size_t)(d * L) * sizeof(R));
41  alpha = (C*) NFFT(malloc)((size_t)(L) * sizeof(C));
42 
44  k = 0;
45  while (k < L)
46  {
47  R r_max = K(1.0);
48  R r2 = K(0.0);
49 
50  for (j = 0; j < d; j++)
51  x[k * d + j] = K(2.0) * r_max * NFFT(drand48)() - r_max;
52 
53  for (j = 0; j < d; j++)
54  r2 += x[k * d + j] * x[k * d + j];
55 
56  if (r2 >= r_max * r_max)
57  continue;
58 
59  k++;
60  }
61 
62  NFFT(vrand_unit_complex)(alpha, L);
63 
65  k = 0;
66  while (k < M)
67  {
68  R r_max = K(1.0);
69  R r2 = K(0.0);
70 
71  for (j = 0; j < d; j++)
72  y[k * d + j] = K(2.0) * r_max * NFFT(drand48)() - r_max;
73 
74  for (j = 0; j < d; j++)
75  r2 += y[k * d + j] * y[k * d + j];
76 
77  if (r2 >= r_max * r_max)
78  continue;
79 
80  k++;
81  }
82 
83  printf("%d %d %d\n", d, L, M);
84 
85  for (j = 0; j < L; j++)
86  {
87  for (t = 0; t < d; t++)
88  printf("%.16" __FES__ " ", x[d * j + t]);
89  printf("\n");
90  }
91 
92  for (j = 0; j < L; j++)
93  printf("%.16" __FES__ " %.16" __FES__ "\n", CREAL(alpha[j]), CIMAG(alpha[j]));
94 
95  for (j = 0; j < M; j++)
96  {
97  for (t = 0; t < d; t++)
98  printf("%.16" __FES__ " ", y[d * j + t]);
99  printf("\n");
100  }
101 
102  NFFT(free)(x);
103  NFFT(free)(y);
104  NFFT(free)(alpha);
105 }
106 
107 int main(int argc, char **argv)
108 {
109  int d;
110  int L;
111  int M;
112 
113  if (argc < 4)
114  {
115  fprintf(stderr, "usage: d L M\n");
116  return -1;
117  }
118 
119  d = atoi(argv[1]);
120  L = atoi(argv[2]);
121  M = atoi(argv[3]);
122 
123  fprintf(stderr, "d=%d, L=%d, M=%d\n", d, L, M);
124 
125  fastsum_benchomp_createdataset(d, L, M);
126 
127  return 0;
128 }
129