NFFT  3.3.1
simple_test.c
00001 /*
00002  * Copyright (c) 2002, 2016 Jens Keiner, Stefan Kunis, Daniel Potts
00003  *
00004  * This program is free software; you can redistribute it and/or modify it under
00005  * the terms of the GNU General Public License as published by the Free Software
00006  * Foundation; either version 2 of the License, or (at your option) any later
00007  * version.
00008  *
00009  * This program is distributed in the hope that it will be useful, but WITHOUT
00010  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00011  * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
00012  * details.
00013  *
00014  * You should have received a copy of the GNU General Public License along with
00015  * this program; if not, write to the Free Software Foundation, Inc., 51
00016  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00017  */
00018 
00019 #include <stdlib.h>
00020 #include <math.h>
00021 #include <complex.h>
00022 
00023 #include "nfft3.h"
00024 
00026 #define CSWAP(x,y) {double _Complex * NFFT_SWAP_temp__; \
00027   NFFT_SWAP_temp__=(x); (x)=(y); (y)=NFFT_SWAP_temp__;}
00028 
00029 void simple_test_nnfft_1d(void)
00030 {
00031   int j,k;                              
00032   nnfft_plan my_plan;                    
00034   int N[1];
00035   N[0]=10;
00036 
00038   nnfft_init(&my_plan, 1, 3, 19, N);
00039 
00041   for(j=0;j<my_plan.M_total;j++)
00042   {
00043     my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5;
00044   }
00046   for(j=0;j<my_plan.N_total;j++)
00047   {
00048     my_plan.v[j]=((double)rand())/((double)RAND_MAX)-0.5;
00049   }
00050 
00052 /*  if(my_plan.nnfft_flags & PRE_PSI)
00053         nnfft_precompute_psi(&my_plan);
00054 
00055   if(my_plan.nnfft_flags & PRE_FULL_PSI)
00056   nnfft_precompute_full_psi(&my_plan);
00057   if(my_plan.nnfft_flags & PRE_LIN_PSI)
00058    nnfft_precompute_lin_psi(&my_plan);
00060 /*  if(my_plan.nnfft_flags & PRE_PHI_HUT)
00061     nnfft_precompute_phi_hut(&my_plan);
00062 */
00063 
00064 nnfft_precompute_one_psi(&my_plan);
00065 
00066 
00068   for(k=0;k<my_plan.N_total;k++)
00069     my_plan.f_hat[k] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
00070 
00071   nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,"given Fourier coefficients, vector f_hat");
00072 
00074   nnfft_trafo_direct(&my_plan);
00075   nfft_vpr_complex(my_plan.f,my_plan.M_total,"nndft, vector f");
00076 
00078   nnfft_trafo(&my_plan);
00079   nfft_vpr_complex(my_plan.f,my_plan.M_total,"nnfft, vector f");
00080 
00082   nnfft_finalize(&my_plan);
00083 }
00084 
00085 static void simple_test_adjoint_nnfft_1d(void)
00086 {
00087   int j;                                 
00088   nnfft_plan my_plan;                    
00090   int N[1];
00091   N[0]=12;
00092 
00094   nnfft_init(&my_plan, 1, 20, 33, N);
00095 
00097   for(j=0;j<my_plan.M_total;j++)
00098   {
00099     my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5;
00100   }
00102   for(j=0;j<my_plan.N_total;j++)
00103   {
00104     my_plan.v[j]=((double)rand())/((double)RAND_MAX)-0.5;
00105   }
00106 
00108   if(my_plan.nnfft_flags & PRE_PSI)
00109     nnfft_precompute_psi(&my_plan);
00110 
00111   if(my_plan.nnfft_flags & PRE_FULL_PSI)
00112     nnfft_precompute_full_psi(&my_plan);
00113 
00114   if(my_plan.nnfft_flags & PRE_LIN_PSI)
00115     nnfft_precompute_lin_psi(&my_plan);
00116 
00118   if(my_plan.nnfft_flags & PRE_PHI_HUT)
00119     nnfft_precompute_phi_hut(&my_plan);
00120 
00122   for(j=0;j<my_plan.M_total;j++)
00123     my_plan.f[j] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
00124 
00125   nfft_vpr_complex(my_plan.f,my_plan.M_total,"given Samples, vector f");
00126 
00128   nnfft_adjoint_direct(&my_plan);
00129   nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,"adjoint nndft, vector f_hat");
00130 
00132   nnfft_adjoint(&my_plan);
00133   nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,"adjoint nnfft, vector f_hat");
00134 
00136   nnfft_finalize(&my_plan);
00137 }
00138 
00139 static void simple_test_nnfft_2d(void)
00140 {
00141   int j,k;                              
00142   nnfft_plan my_plan;                    
00144   int N[2];
00145   N[0]=12;
00146   N[1]=14;
00147 
00149   nnfft_init(&my_plan, 2,12*14,19, N);
00150 
00152   for(j=0;j<my_plan.M_total;j++)
00153   {
00154     my_plan.x[2*j]=((double)rand())/((double)RAND_MAX)-0.5;
00155     my_plan.x[2*j+1]=((double)rand())/((double)RAND_MAX)-0.5;
00156   }
00157 
00159   for(j=0;j<my_plan.N_total;j++)
00160   {
00161     my_plan.v[2*j]=((double)rand())/((double)RAND_MAX)-0.5;
00162     my_plan.v[2*j+1]=((double)rand())/((double)RAND_MAX)-0.5;
00163   }
00164 
00166   if(my_plan.nnfft_flags & PRE_PSI)
00167     nnfft_precompute_psi(&my_plan);
00168 
00169   if(my_plan.nnfft_flags & PRE_FULL_PSI)
00170     nnfft_precompute_full_psi(&my_plan);
00171 
00172   if(my_plan.nnfft_flags & PRE_LIN_PSI)
00173     nnfft_precompute_lin_psi(&my_plan);
00174 
00176   if(my_plan.nnfft_flags & PRE_PHI_HUT)
00177     nnfft_precompute_phi_hut(&my_plan);
00178 
00180   for(k=0;k<my_plan.N_total;k++)
00181     my_plan.f_hat[k] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
00182 
00183   nfft_vpr_complex(my_plan.f_hat,12,
00184         "given Fourier coefficients, vector f_hat (first 12 entries)");
00185 
00187   nnfft_trafo_direct(&my_plan);
00188   nfft_vpr_complex(my_plan.f,my_plan.M_total,"ndft, vector f");
00189 
00191   nnfft_trafo(&my_plan);
00192   nfft_vpr_complex(my_plan.f,my_plan.M_total,"nfft, vector f");
00193 
00195   nnfft_finalize(&my_plan);
00196 }
00197 
00198 static void simple_test_innfft_1d(void)
00199 {
00200   int j,k,l,N=8;                        
00201   nnfft_plan my_plan;                   
00202   solver_plan_complex my_iplan;         
00205   nnfft_init(&my_plan,1 ,8 ,8 ,&N);
00206 
00208   solver_init_advanced_complex(&my_iplan,(nfft_mv_plan_complex*)(&my_plan),CGNR);
00209 
00211   for(j=0;j<my_plan.M_total;j++)
00212     my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5;
00213 
00215   for(k=0;k<my_plan.N_total;k++)
00216     my_plan.v[k]=((double)rand())/((double)RAND_MAX)-0.5;
00217 
00219   if(my_plan.nnfft_flags & PRE_PSI)
00220     nnfft_precompute_psi(&my_plan);
00221 
00222   if(my_plan.nnfft_flags & PRE_FULL_PSI)
00223       nnfft_precompute_full_psi(&my_plan);
00224 
00226   if(my_plan.nnfft_flags & PRE_PHI_HUT)
00227     nnfft_precompute_phi_hut(&my_plan);
00228 
00230   for(j=0;j<my_plan.M_total;j++)
00231     my_iplan.y[j] = ((double)rand())/((double)RAND_MAX);
00232 
00233   nfft_vpr_complex(my_iplan.y,my_plan.M_total,"given data, vector given_f");
00234 
00236   for(k=0;k<my_plan.N_total;k++)
00237     my_iplan.f_hat_iter[k] = 0.0;
00238 
00239   nfft_vpr_complex(my_iplan.f_hat_iter,my_plan.N_total,
00240         "approximate solution, vector f_hat_iter");
00241 
00243   solver_before_loop_complex(&my_iplan);
00244 
00245   for(l=0;l<8;l++)
00246   {
00247     printf("iteration l=%d\n",l);
00248     solver_loop_one_step_complex(&my_iplan);
00249     nfft_vpr_complex(my_iplan.f_hat_iter,my_plan.N_total,
00250           "approximate solution, vector f_hat_iter");
00251 
00252     CSWAP(my_iplan.f_hat_iter,my_plan.f_hat);
00253     nnfft_trafo(&my_plan);
00254     nfft_vpr_complex(my_plan.f,my_plan.M_total,"fitting the data, vector f");
00255     CSWAP(my_iplan.f_hat_iter,my_plan.f_hat);
00256   }
00257 
00258   solver_finalize_complex(&my_iplan);
00259   nnfft_finalize(&my_plan);
00260 }
00261 
00262 static void measure_time_nnfft_1d(void)
00263 {
00264   int j,k;                              
00265   nnfft_plan my_plan;                    
00266   int my_N,N=4;
00267   double t;
00268   double t0, t1;
00269 
00270   for(my_N=16; my_N<=16384; my_N*=2)
00271   {
00272     nnfft_init(&my_plan,1,my_N,my_N,&N);
00273 
00274     for(j=0;j<my_plan.M_total;j++)
00275       my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5;
00276 
00277     for(j=0;j<my_plan.N_total;j++)
00278       my_plan.v[j]=((double)rand())/((double)RAND_MAX)-0.5;
00279 
00280     if(my_plan.nnfft_flags & PRE_PSI)
00281       nnfft_precompute_psi(&my_plan);
00282 
00283     if(my_plan.nnfft_flags & PRE_FULL_PSI)
00284         nnfft_precompute_full_psi(&my_plan);
00285 
00286     if(my_plan.nnfft_flags & PRE_PHI_HUT)
00287       nnfft_precompute_phi_hut(&my_plan);
00288 
00289     for(k=0;k<my_plan.N_total;k++)
00290       my_plan.f_hat[k] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
00291 
00292     t0 = nfft_clock_gettime_seconds();
00293     nnfft_trafo_direct(&my_plan);
00294     t1 = nfft_clock_gettime_seconds();
00295     t = t1 - t0;
00296     printf("t_nndft=%e,\t",t);
00297 
00298     t0 = nfft_clock_gettime_seconds();
00299     nnfft_trafo(&my_plan);
00300     t1 = nfft_clock_gettime_seconds();
00301     t = t1 - t0;
00302     printf("t_nnfft=%e\t",t);
00303 
00304     printf("(N=M=%d)\n",my_N);
00305 
00306     nnfft_finalize(&my_plan);
00307   }
00308 }
00309 
00310 int main(void)
00311 {
00312   system("clear");
00313   printf("1) computing a one dimensional nndft, nnfft\n\n");
00314   simple_test_nnfft_1d();
00315 
00316   /*getc(stdin);
00317 
00318   system("clear");
00319   printf("1a) computing a one dimensional adjoint nndft, nnfft\n\n");
00320   simple_test_adjoint_nnfft_1d();
00321 
00322   getc(stdin);
00323 
00324   system("clear");
00325   printf("2) computing a two dimensional nndft, nnfft\n\n");
00326   simple_test_nnfft_2d();
00327 
00328   getc(stdin);
00329 
00330   system("clear");
00331   printf("3) computing a one dimensional innfft\n\n");
00332   simple_test_innfft_1d();
00333 
00334   getc(stdin);
00335 
00336   system("clear");
00337   printf("4) computing times for one dimensional nnfft\n\n");
00338   measure_time_nnfft_1d();
00339 */
00340   return 1;
00341 }