NFFT Logo 3.2.2
solver/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 <stdio.h>
00023 #include <math.h>
00024 #include <string.h>
00025 #include <stdlib.h>
00026 #ifdef HAVE_COMPLEX_H
00027 #include <complex.h>
00028 #endif
00029 
00030 #include "nfft3util.h"
00031 #include "nfft3.h"
00032 
00033 /* void simple_test_infft_1d(int N, int M, int iter) */
00034 /* { */
00035 /*   int k,l;                            /\**< index for nodes, freqencies,iter*\/ */
00036 /*   nfft_plan p;                          /\**< plan for the nfft               *\/ */
00037 /*   infft_plan ip;                        /\**< plan for the inverse nfft       *\/ */
00038 
00039 /*   /\** initialise an one dimensional plan *\/ */
00040 /*   nfft_init_1d(&p, N, M); */
00041 
00042 /*   /\** init pseudo random nodes *\/ */
00043 /*   nfft_vrand_shifted_unit_double(p.x,p.M_total); */
00044 
00045 /*   /\** precompute psi, the entries of the matrix B *\/ */
00046 /*   if(p.nfft_flags & PRE_ONE_PSI) */
00047 /*     nfft_precompute_one_psi(&p); */
00048 
00049 /*   /\** initialise inverse plan *\/ */
00050 /*   infft_init(&ip,&p); */
00051 
00052 /*   /\** init pseudo random samples and show them *\/ */
00053 /*   nfft_vrand_unit_complex(ip.y,p.M_total); */
00054 /*   nfft_vpr_complex(ip.y,p.M_total,"Given data, vector y"); */
00055 
00056 /*   /\** initialise some guess f_hat_0 and solve *\/ */
00057 /*   for(k=0;k<p.N_total;k++) */
00058 /*       ip.f_hat_iter[k]=0; */
00059 
00060 /*   nfft_vpr_complex(ip.f_hat_iter,p.N_total,"Initial guess, vector f_hat_iter"); */
00061 
00062 /*   NFFT_SWAP_complex(ip.f_hat_iter,p.f_hat); */
00063 /*   nfft_trafo(&p); */
00064 /*   nfft_vpr_complex(p.f,p.M_total,"Data fit, vector f"); */
00065 /*   NFFT_SWAP_complex(ip.f_hat_iter,p.f_hat); */
00066 
00067 /*   infft_before_loop(&ip); */
00068 /*   printf("\n Residual r=%e\n",ip.dot_r_iter); */
00069 
00070 /*   for(l=0;l<iter;l++) */
00071 /*     { */
00072 /*       printf("\n********** Iteration l=%d **********\n",l); */
00073 /*       infft_loop_one_step(&ip); */
00074 /*       nfft_vpr_complex(ip.f_hat_iter,p.N_total, */
00075 /*      "Approximate solution, vector f_hat_iter"); */
00076 
00077 /*       NFFT_SWAP_complex(ip.f_hat_iter,p.f_hat); */
00078 /*       nfft_trafo(&p); */
00079 /*       nfft_vpr_complex(p.f,p.M_total,"Data fit, vector f"); */
00080 /*       NFFT_SWAP_complex(ip.f_hat_iter,p.f_hat); */
00081 
00082 /*       printf("\n Residual r=%e\n",ip.dot_r_iter); */
00083 /*     } */
00084 
00085 /*   infft_finalize(&ip); */
00086 /*   nfft_finalize(&p); */
00087 /* } */
00088 
00090 static void simple_test_solver_nfft_1d(int N, int M, int iter)
00091 {
00092   int k,l;                            
00093   nfft_plan p;                          
00094   solver_plan_complex ip;                        
00097   nfft_init_1d(&p, N, M);
00098 
00100   nfft_vrand_shifted_unit_double(p.x,p.M_total);
00101 
00103   if(p.nfft_flags & PRE_ONE_PSI)
00104     nfft_precompute_one_psi(&p);
00105 
00107   solver_init_complex(&ip,(nfft_mv_plan_complex*)(&p));
00108 
00110   nfft_vrand_unit_complex(ip.y,p.M_total);
00111   nfft_vpr_complex(ip.y,p.M_total,"Given data, vector y");
00112 
00114   for(k=0;k<p.N_total;k++)
00115       ip.f_hat_iter[k]=0;
00116 
00117   nfft_vpr_complex(ip.f_hat_iter,p.N_total,"Initial guess, vector f_hat_iter");
00118 
00119   NFFT_SWAP_complex(ip.f_hat_iter,p.f_hat);
00120   nfft_trafo(&p);
00121   nfft_vpr_complex(p.f,p.M_total,"Data fit, vector f");
00122   NFFT_SWAP_complex(ip.f_hat_iter,p.f_hat);
00123 
00124   solver_before_loop_complex(&ip);
00125   printf("\n Residual r=%e\n",ip.dot_r_iter);
00126 
00127   for(l=0;l<iter;l++)
00128     {
00129       printf("\n********** Iteration l=%d **********\n",l);
00130       solver_loop_one_step_complex(&ip);
00131       nfft_vpr_complex(ip.f_hat_iter,p.N_total,
00132       "Approximate solution, vector f_hat_iter");
00133 
00134       NFFT_SWAP_complex(ip.f_hat_iter,p.f_hat);
00135       nfft_trafo(&p);
00136       nfft_vpr_complex(p.f,p.M_total,"Data fit, vector f");
00137       NFFT_SWAP_complex(ip.f_hat_iter,p.f_hat);
00138 
00139       printf("\n Residual r=%e\n",ip.dot_r_iter);
00140     }
00141 
00142   solver_finalize_complex(&ip);
00143   nfft_finalize(&p);
00144 }
00145 
00147 int main(void)
00148 {
00149   printf("\n Computing a one dimensional inverse nfft\n");
00150 
00151   simple_test_solver_nfft_1d(8,4,5);
00152 
00153   return 1;
00154 }
00155 /* \} */

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409