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 nfsft_benchomp_createdataset(unsigned int trafo_adjoint, int N, int M)
00034 {
00035 int t, j, k, n;
00036 R *x;
00037 C *f, *f_hat;
00038 int N_total = (2*N+2) * (2*N+2);
00039 nfsft_plan ptemp;
00040
00041 nfsft_init_guru(&ptemp, N, M, NFSFT_MALLOC_X | NFSFT_MALLOC_F |
00042 NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_PRESERVE_F_HAT,
00043 PRE_PHI_HUT | PRE_PSI | FFTW_INIT | FFT_OUT_OF_PLACE, 6);
00044
00045 x = (R*) nfft_malloc(2*M*sizeof(R));
00046 f = (C*) nfft_malloc(M*sizeof(C));
00047 f_hat = (C*) nfft_malloc(N_total*sizeof(C));
00048
00049
00050 for (j = 0; j < M; j++)
00051 {
00052 x[2*j]= nfft_drand48() - K(0.5);
00053 x[2*j+1]= K(0.5) * nfft_drand48();
00054 }
00055
00056 if (trafo_adjoint==0)
00057 {
00058 for (k = 0; k <= N; k++)
00059 for (n = -k; n <= k; n++)
00060 nfft_vrand_unit_complex(f_hat+NFSFT_INDEX(k,n,&ptemp),1);
00061 }
00062 else
00063 {
00064 nfft_vrand_unit_complex(f,M);
00065 }
00066
00067 printf("%d %d %d\n", trafo_adjoint, N, M);
00068
00069 for (j=0; j < M; j++)
00070 {
00071 for (t=0; t < 2; t++)
00072 printf("%.16e ", x[2*j+t]);
00073 printf("\n");
00074 }
00075
00076 if (trafo_adjoint==0)
00077 {
00078 for (k = 0; k <= N; k++)
00079 for (n = -k; n <= k; n++)
00080 printf("%.16e %.16e\n", creal(f_hat[NFSFT_INDEX(k,n,&ptemp)]), cimag(f_hat[NFSFT_INDEX(k,n,&ptemp)]));
00081 }
00082 else
00083 {
00084 for (j=0; j < M; j++)
00085 printf("%.16e %.16e\n", creal(f[j]), cimag(f[j]));
00086 }
00087
00088 nfft_free(x);
00089 nfft_free(f);
00090 nfft_free(f_hat);
00091 }
00092
00093 int main(int argc, char **argv)
00094 {
00095 int trafo_adjoint;
00096 int N;
00097 int M;
00098
00099 if (argc < 4) {
00100 fprintf(stderr, "usage: tr_adj N M\n");
00101 return -1;
00102 }
00103
00104 trafo_adjoint = atoi(argv[1]);
00105 if (trafo_adjoint < 0 && trafo_adjoint > 1)
00106 trafo_adjoint = 1;
00107
00108 N = atoi(argv[2]);
00109 M = atoi(argv[3]);
00110 fprintf(stderr, "tr_adj=%d, N=%d, M=%d\n", trafo_adjoint, N, M);
00111
00112 nfsft_benchomp_createdataset(trafo_adjoint, N, M);
00113
00114 return 0;
00115 }