NFFT Logo 3.2.2
nnfft/simple_test.c
00001 /*
00002  * Copyright (c) 2002, 2012 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 /* $Id: simple_test.c 3896 2012-10-10 12:19:26Z tovo $ */
00020 #include "config.h"
00021 
00022 #include <stdlib.h>
00023 #include <math.h>
00024 #ifdef HAVE_COMPLEX_H
00025 #include <complex.h>
00026 #endif
00027 
00028 #include "nfft3util.h"
00029 #include "nfft3.h"
00030 #include "infft.h"
00031 
00032 void simple_test_nnfft_1d(void)
00033 {
00034   int j,k;                              
00035   nnfft_plan my_plan;                    
00037   int N[1];
00038   N[0]=12;
00039 
00041   nnfft_init(&my_plan, 1, 3, 19, N);
00042 
00044   for(j=0;j<my_plan.M_total;j++)
00045   {
00046     my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5;
00047   }
00049   for(j=0;j<my_plan.N_total;j++)
00050   {
00051     my_plan.v[j]=((double)rand())/((double)RAND_MAX)-0.5;
00052   }
00053 
00055   if(my_plan.nnfft_flags & PRE_PSI)
00056     nnfft_precompute_psi(&my_plan);
00057 
00058   if(my_plan.nnfft_flags & PRE_FULL_PSI)
00059     nnfft_precompute_full_psi(&my_plan);
00060 
00061   if(my_plan.nnfft_flags & PRE_LIN_PSI)
00062     nnfft_precompute_lin_psi(&my_plan);
00063 
00065   if(my_plan.nnfft_flags & PRE_PHI_HUT)
00066     nnfft_precompute_phi_hut(&my_plan);
00067 
00069   for(k=0;k<my_plan.N_total;k++)
00070     my_plan.f_hat[k] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
00071 
00072   nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,"given Fourier coefficients, vector f_hat");
00073 
00075   nnfft_trafo_direct(&my_plan);
00076   nfft_vpr_complex(my_plan.f,my_plan.M_total,"nndft, vector f");
00077 
00079   nnfft_trafo(&my_plan);
00080   nfft_vpr_complex(my_plan.f,my_plan.M_total,"nnfft, vector f");
00081 
00083   nnfft_finalize(&my_plan);
00084 }
00085 
00086 static void simple_test_adjoint_nnfft_1d(void)
00087 {
00088   int j;                                 
00089   nnfft_plan my_plan;                    
00091   int N[1];
00092   N[0]=12;
00093 
00095   nnfft_init(&my_plan, 1, 20, 33, N);
00096 
00098   for(j=0;j<my_plan.M_total;j++)
00099   {
00100     my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5;
00101   }
00103   for(j=0;j<my_plan.N_total;j++)
00104   {
00105     my_plan.v[j]=((double)rand())/((double)RAND_MAX)-0.5;
00106   }
00107 
00109   if(my_plan.nnfft_flags & PRE_PSI)
00110     nnfft_precompute_psi(&my_plan);
00111 
00112   if(my_plan.nnfft_flags & PRE_FULL_PSI)
00113     nnfft_precompute_full_psi(&my_plan);
00114 
00115   if(my_plan.nnfft_flags & PRE_LIN_PSI)
00116     nnfft_precompute_lin_psi(&my_plan);
00117 
00119   if(my_plan.nnfft_flags & PRE_PHI_HUT)
00120     nnfft_precompute_phi_hut(&my_plan);
00121 
00123   for(j=0;j<my_plan.M_total;j++)
00124     my_plan.f[j] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
00125 
00126   nfft_vpr_complex(my_plan.f,my_plan.M_total,"given Samples, vector f");
00127 
00129   nnfft_adjoint_direct(&my_plan);
00130   nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,"adjoint nndft, vector f_hat");
00131 
00133   nnfft_adjoint(&my_plan);
00134   nfft_vpr_complex(my_plan.f_hat,my_plan.N_total,"adjoint nnfft, vector f_hat");
00135 
00137   nnfft_finalize(&my_plan);
00138 }
00139 
00140 static void simple_test_nnfft_2d(void)
00141 {
00142   int j,k;                              
00143   nnfft_plan my_plan;                    
00145   int N[2];
00146   N[0]=12;
00147   N[1]=14;
00148 
00150   nnfft_init(&my_plan, 2,12*14,19, N);
00151 
00153   for(j=0;j<my_plan.M_total;j++)
00154   {
00155     my_plan.x[2*j]=((double)rand())/((double)RAND_MAX)-0.5;
00156     my_plan.x[2*j+1]=((double)rand())/((double)RAND_MAX)-0.5;
00157   }
00158 
00160   for(j=0;j<my_plan.N_total;j++)
00161   {
00162     my_plan.v[2*j]=((double)rand())/((double)RAND_MAX)-0.5;
00163     my_plan.v[2*j+1]=((double)rand())/((double)RAND_MAX)-0.5;
00164   }
00165 
00167   if(my_plan.nnfft_flags & PRE_PSI)
00168     nnfft_precompute_psi(&my_plan);
00169 
00170   if(my_plan.nnfft_flags & PRE_FULL_PSI)
00171     nnfft_precompute_full_psi(&my_plan);
00172 
00173   if(my_plan.nnfft_flags & PRE_LIN_PSI)
00174     nnfft_precompute_lin_psi(&my_plan);
00175 
00177   if(my_plan.nnfft_flags & PRE_PHI_HUT)
00178     nnfft_precompute_phi_hut(&my_plan);
00179 
00181   for(k=0;k<my_plan.N_total;k++)
00182     my_plan.f_hat[k] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
00183 
00184   nfft_vpr_complex(my_plan.f_hat,12,
00185         "given Fourier coefficients, vector f_hat (first 12 entries)");
00186 
00188   nnfft_trafo_direct(&my_plan);
00189   nfft_vpr_complex(my_plan.f,my_plan.M_total,"ndft, vector f");
00190 
00192   nnfft_trafo(&my_plan);
00193   nfft_vpr_complex(my_plan.f,my_plan.M_total,"nfft, vector f");
00194 
00196   nnfft_finalize(&my_plan);
00197 }
00198 
00199 static void simple_test_innfft_1d(void)
00200 {
00201   int j,k,l,N=8;                        
00202   nnfft_plan my_plan;                   
00203   solver_plan_complex my_iplan;         
00206   nnfft_init(&my_plan,1 ,8 ,8 ,&N);
00207 
00209   solver_init_advanced_complex(&my_iplan,(nfft_mv_plan_complex*)(&my_plan),CGNR);
00210 
00212   for(j=0;j<my_plan.M_total;j++)
00213     my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5;
00214 
00216   for(k=0;k<my_plan.N_total;k++)
00217     my_plan.v[k]=((double)rand())/((double)RAND_MAX)-0.5;
00218 
00220   if(my_plan.nnfft_flags & PRE_PSI)
00221     nnfft_precompute_psi(&my_plan);
00222 
00223   if(my_plan.nnfft_flags & PRE_FULL_PSI)
00224       nnfft_precompute_full_psi(&my_plan);
00225 
00227   if(my_plan.nnfft_flags & PRE_PHI_HUT)
00228     nnfft_precompute_phi_hut(&my_plan);
00229 
00231   for(j=0;j<my_plan.M_total;j++)
00232     my_iplan.y[j] = ((double)rand())/((double)RAND_MAX);
00233 
00234   nfft_vpr_complex(my_iplan.y,my_plan.M_total,"given data, vector given_f");
00235 
00237   for(k=0;k<my_plan.N_total;k++)
00238     my_iplan.f_hat_iter[k] = 0.0;
00239 
00240   nfft_vpr_complex(my_iplan.f_hat_iter,my_plan.N_total,
00241         "approximate solution, vector f_hat_iter");
00242 
00244   solver_before_loop_complex(&my_iplan);
00245 
00246   for(l=0;l<8;l++)
00247   {
00248     printf("iteration l=%d\n",l);
00249     solver_loop_one_step_complex(&my_iplan);
00250     nfft_vpr_complex(my_iplan.f_hat_iter,my_plan.N_total,
00251           "approximate solution, vector f_hat_iter");
00252 
00253     NFFT_SWAP_complex(my_iplan.f_hat_iter,my_plan.f_hat);
00254     nnfft_trafo(&my_plan);
00255     nfft_vpr_complex(my_plan.f,my_plan.M_total,"fitting the data, vector f");
00256     NFFT_SWAP_complex(my_iplan.f_hat_iter,my_plan.f_hat);
00257   }
00258 
00259   solver_finalize_complex(&my_iplan);
00260   nnfft_finalize(&my_plan);
00261 }
00262 
00263 static void measure_time_nnfft_1d(void)
00264 {
00265   int j,k;                              
00266   nnfft_plan my_plan;                    
00267   int my_N,N=4;
00268   double t;
00269   ticks t0, t1;
00270 
00271   for(my_N=16; my_N<=16384; my_N*=2)
00272   {
00273     nnfft_init(&my_plan,1,my_N,my_N,&N);
00274 
00275     for(j=0;j<my_plan.M_total;j++)
00276       my_plan.x[j]=((double)rand())/((double)RAND_MAX)-0.5;
00277 
00278     for(j=0;j<my_plan.N_total;j++)
00279       my_plan.v[j]=((double)rand())/((double)RAND_MAX)-0.5;
00280 
00281     if(my_plan.nnfft_flags & PRE_PSI)
00282       nnfft_precompute_psi(&my_plan);
00283 
00284     if(my_plan.nnfft_flags & PRE_FULL_PSI)
00285         nnfft_precompute_full_psi(&my_plan);
00286 
00287     if(my_plan.nnfft_flags & PRE_PHI_HUT)
00288       nnfft_precompute_phi_hut(&my_plan);
00289 
00290     for(k=0;k<my_plan.N_total;k++)
00291       my_plan.f_hat[k] = ((double)rand())/((double)RAND_MAX) + _Complex_I*((double)rand())/((double)RAND_MAX);
00292 
00293     t0 = getticks();
00294     nnfft_trafo_direct(&my_plan);
00295     t1 = getticks();
00296     t = nfft_elapsed_seconds(t1,t0);
00297     printf("t_nndft=%e,\t",t);
00298 
00299     t0 = getticks();
00300     nnfft_trafo(&my_plan);
00301     t1 = getticks();
00302     t = nfft_elapsed_seconds(t1,t0);
00303     printf("t_nnfft=%e\t",t);
00304 
00305     printf("(N=M=%d)\n",my_N);
00306 
00307     nnfft_finalize(&my_plan);
00308   }
00309 }
00310 
00311 int main(void)
00312 {
00313   system("clear");
00314   printf("1) computing a one dimensional nndft, nnfft\n\n");
00315   simple_test_nnfft_1d();
00316 
00317   getc(stdin);
00318 
00319   system("clear");
00320   printf("1a) computing a one dimensional adjoint nndft, nnfft\n\n");
00321   simple_test_adjoint_nnfft_1d();
00322 
00323   getc(stdin);
00324 
00325   system("clear");
00326   printf("2) computing a two dimensional nndft, nnfft\n\n");
00327   simple_test_nnfft_2d();
00328 
00329   getc(stdin);
00330 
00331   system("clear");
00332   printf("3) computing a one dimensional innfft\n\n");
00333   simple_test_innfft_1d();
00334 
00335   getc(stdin);
00336 
00337   system("clear");
00338   printf("4) computing times for one dimensional nnfft\n\n");
00339   measure_time_nnfft_1d();
00340 
00341   return 1;
00342 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409