00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "config.h"
00022
00023
00024 #include <stdio.h>
00025 #include <math.h>
00026 #include <string.h>
00027 #include <stdlib.h>
00028
00029 #ifdef HAVE_COMPLEX_H
00030 #include <complex.h>
00031 #endif
00032 #include <omp.h>
00033
00034 #include "nfft3.h"
00035 #include "nfft3util.h"
00036 #include "infft.h"
00037
00038 static void simple_test_nfsft(void)
00039 {
00040 const int N = 4;
00041 const int M = 8;
00042 nfsft_plan plan;
00043 int j, k, n;
00044
00045
00046 nfsft_precompute(N,1000.0,0U,0U);
00047
00048
00049
00050
00051
00052
00053
00054 nfsft_init_guru(&plan, N, M, NFSFT_MALLOC_X | NFSFT_MALLOC_F |
00055 NFSFT_MALLOC_F_HAT | NFSFT_NORMALIZED | NFSFT_PRESERVE_F_HAT,
00056 PRE_PHI_HUT | PRE_PSI | FFTW_INIT | FFT_OUT_OF_PLACE, 6);
00057
00058
00059 for (j = 0; j < plan.M_total; j++)
00060 {
00061 plan.x[2*j]= nfft_drand48() - K(0.5);
00062 plan.x[2*j+1]= K(0.5) * nfft_drand48();
00063 }
00064
00065
00066 nfsft_precompute_x(&plan);
00067
00068
00069 for (k = 0; k <= plan.N; k++)
00070 for (n = -k; n <= k; n++)
00071 plan.f_hat[NFSFT_INDEX(k,n,&plan)] =
00072 nfft_drand48() - K(0.5) + _Complex_I*(nfft_drand48() - K(0.5));
00073
00074
00075 nfsft_trafo_direct(&plan);
00076 printf("Vector f (NDSFT):\n");
00077 for (j = 0; j < plan.M_total; j++)
00078 printf("f[%+2d] = %+5.3" FE " %+5.3" FE "*I\n",j,
00079 creal(plan.f[j]), cimag(plan.f[j]));
00080
00081 printf("\n");
00082
00083
00084 printf("Vector f (NDSFT):\n");
00085 for (j = 0; j < plan.M_total; j++)
00086 printf("f[%+2d] = %+5.3" FE " %+5.3" FE "*I\n",j,
00087 creal(plan.f[j]), cimag(plan.f[j]));
00088
00089 printf("\n");
00090
00091
00092 nfsft_adjoint_direct(&plan);
00093 printf("Vector f_hat (NDSFT):\n");
00094 for (k = 0; k <= plan.N; k++)
00095 for (n = -k; n <= k; n++)
00096 fprintf(stdout,"f_hat[%+2d,%+2d] = %+5.3" FE " %+5.3" FE "*I\n",k,n,
00097 creal(plan.f_hat[NFSFT_INDEX(k,n,&plan)]),
00098 cimag(plan.f_hat[NFSFT_INDEX(k,n,&plan)]));
00099
00100 printf("\n");
00101
00102
00103 nfsft_adjoint(&plan);
00104 printf("Vector f_hat (NFSFT):\n");
00105 for (k = 0; k <= plan.N; k++)
00106 {
00107 for (n = -k; n <= k; n++)
00108 {
00109 fprintf(stdout,"f_hat[%+2d,%+2d] = %+5.3" FE " %+5.3" FE "*I\n",k,n,
00110 creal(plan.f_hat[NFSFT_INDEX(k,n,&plan)]),
00111 cimag(plan.f_hat[NFSFT_INDEX(k,n,&plan)]));
00112 }
00113 }
00114
00115
00116 nfsft_finalize(&plan);
00117
00118
00119 nfsft_forget();
00120 }
00121
00122 int main(void)
00123 {
00124 printf("nthreads = %d\n", nfft_get_omp_num_threads());
00125
00126
00127 fftw_init_threads();
00128
00129 printf("Computing an NDSFT, an NFSFT, an adjoint NDSFT, and an adjoint NFSFT"
00130 "...\n\n");
00131 simple_test_nfsft();
00132 return EXIT_SUCCESS;
00133 }