00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include <stdio.h>
00022 #include <math.h>
00023 #include <string.h>
00024 #include <stdlib.h>
00025 #include <complex.h>
00026
00027 #include "config.h"
00028
00029 #include "nfft3util.h"
00030 #include "nfft3.h"
00031 #include "infft.h"
00032
00033 void nfft_benchomp_createdataset(unsigned int d, unsigned int trafo_adjoint, int *N, int M, double sigma)
00034 {
00035 int n[d];
00036 int t, j;
00037 R *x;
00038 C *f, *f_hat;
00039 int N_total = 1;
00040
00041 for (t = 0; t < d; t++)
00042 N_total *= N[t];
00043
00044 x = (R*) nfft_malloc(d*M*sizeof(R));
00045 f = (C*) nfft_malloc(M*sizeof(C));
00046 f_hat = (C*) nfft_malloc(N_total*sizeof(C));
00047
00048 for (t=0; t<d; t++)
00049 n[t] = sigma*nfft_next_power_of_2(N[t]);
00050
00052 nfft_vrand_shifted_unit_double(x,d*M);
00053
00054 if (trafo_adjoint==0)
00055 {
00056 nfft_vrand_unit_complex(f_hat,N_total);
00057 }
00058 else
00059 {
00060 nfft_vrand_unit_complex(f,M);
00061 }
00062
00063 printf("%d %d ", d, trafo_adjoint);
00064
00065 for (t=0; t<d; t++)
00066 printf("%d ", N[t]);
00067
00068 for (t=0; t<d; t++)
00069 printf("%d ", n[t]);
00070
00071 printf("%d\n", M);
00072
00073 for (j=0; j < M; j++)
00074 {
00075 for (t=0; t < d; t++)
00076 printf("%.16e ", x[d*j+t]);
00077 printf("\n");
00078 }
00079
00080 if (trafo_adjoint==0)
00081 {
00082 for (j=0; j < N_total; j++)
00083 printf("%.16e %.16e\n", creal(f_hat[j]), cimag(f_hat[j]));
00084 }
00085 else
00086 {
00087 for (j=0; j < M; j++)
00088 printf("%.16e %.16e\n", creal(f[j]), cimag(f[j]));
00089 }
00090
00091 nfft_free(x);
00092 nfft_free(f);
00093 nfft_free(f_hat);
00094 }
00095
00096 int main(int argc, char **argv)
00097 {
00098 int d;
00099 int *N;
00100 int M;
00101 int t;
00102 int trafo_adjoint;
00103 double sigma;
00104
00105 if (argc < 6) {
00106 fprintf(stderr, "usage: d tr_adj N_1 ... N_d M sigma\n");
00107 return -1;
00108 }
00109
00110 d = atoi(argv[1]);
00111
00112 fprintf(stderr, "d=%d", d);
00113
00114 if (d < 1 || argc < 5+d) {
00115 fprintf(stderr, "usage: d tr_adj N_1 ... N_d M sigma\n");
00116 return -1;
00117 }
00118
00119 N = malloc(d*sizeof(int));
00120
00121 trafo_adjoint = atoi(argv[2]);
00122 if (trafo_adjoint < 0 && trafo_adjoint > 1)
00123 trafo_adjoint = 1;
00124
00125 fprintf(stderr, ", tr_adj=%d, N=", trafo_adjoint);
00126
00127 for (t=0; t<d; t++)
00128 N[t] = atoi(argv[3+t]);
00129
00130 for (t=0; t<d; t++)
00131 fprintf(stderr, "%d ",N[t]);
00132
00133
00134 M = atoi(argv[3+d]);
00135 sigma = atof(argv[4+d]);
00136
00137 fprintf(stderr, ", M=%d, sigma=%.16g\n", M, sigma);
00138
00139 nfft_benchomp_createdataset(d, trafo_adjoint, N, M, sigma);
00140
00141 free(N);
00142
00143 return 0;
00144 }