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 <stdio.h>
00020 #include <math.h>
00021 #include <string.h>
00022 #include <stdlib.h>
00023 #include <complex.h>
00024 
00025 #define NFFT_PRECISION_DOUBLE
00026 
00027 #include "nfft3mp.h"
00028 
00029 /* void simple_test_infft_1d(int N, int M, int iter) */
00030 /* { */
00031 /*   int k,l;                            /\**< index for nodes, freqencies,iter*\/ */
00032 /*   nfft_plan p;                          /\**< plan for the nfft               *\/ */
00033 /*   infft_plan ip;                        /\**< plan for the inverse nfft       *\/ */
00034 
00035 /*   /\** initialise an one dimensional plan *\/ */
00036 /*   nfft_init_1d(&p, N, M); */
00037 
00038 /*   /\** init pseudo random nodes *\/ */
00039 /*   nfft_vrand_shifted_unit_double(p.x,p.M_total); */
00040 
00041 /*   /\** precompute psi, the entries of the matrix B *\/ */
00042 /*   if(p.flags & PRE_ONE_PSI) */
00043 /*     nfft_precompute_one_psi(&p); */
00044 
00045 /*   /\** initialise inverse plan *\/ */
00046 /*   infft_init(&ip,&p); */
00047 
00048 /*   /\** init pseudo random samples and show them *\/ */
00049 /*   nfft_vrand_unit_complex(ip.y,p.M_total); */
00050 /*   nfft_vpr_complex(ip.y,p.M_total,"Given data, vector y"); */
00051 
00052 /*   /\** initialise some guess f_hat_0 and solve *\/ */
00053 /*   for(k=0;k<p.N_total;k++) */
00054 /*       ip.f_hat_iter[k]=0; */
00055 
00056 /*   nfft_vpr_complex(ip.f_hat_iter,p.N_total,"Initial guess, vector f_hat_iter"); */
00057 
00058 /*   CSWAP(ip.f_hat_iter,p.f_hat); */
00059 /*   nfft_trafo(&p); */
00060 /*   nfft_vpr_complex(p.f,p.M_total,"Data fit, vector f"); */
00061 /*   CSWAP(ip.f_hat_iter,p.f_hat); */
00062 
00063 /*   infft_before_loop(&ip); */
00064 /*   printf("\n Residual r=%e\n",ip.dot_r_iter); */
00065 
00066 /*   for(l=0;l<iter;l++) */
00067 /*     { */
00068 /*       printf("\n********** Iteration l=%d **********\n",l); */
00069 /*       infft_loop_one_step(&ip); */
00070 /*       nfft_vpr_complex(ip.f_hat_iter,p.N_total, */
00071 /*      "Approximate solution, vector f_hat_iter"); */
00072 
00073 /*       CSWAP(ip.f_hat_iter,p.f_hat); */
00074 /*       nfft_trafo(&p); */
00075 /*       nfft_vpr_complex(p.f,p.M_total,"Data fit, vector f"); */
00076 /*       CSWAP(ip.f_hat_iter,p.f_hat); */
00077 
00078 /*       printf("\n Residual r=%e\n",ip.dot_r_iter); */
00079 /*     } */
00080 
00081 /*   infft_finalize(&ip); */
00082 /*   nfft_finalize(&p); */
00083 /* } */
00084 
00086 static void simple_test_solver_nfft_1d(int N, int M, int iter)
00087 {
00088   int k, l; 
00089   NFFT(plan) p; 
00090   SOLVER(plan_complex) ip; 
00091   const char *error_str;
00092 
00094   NFFT(init_1d)(&p, N, M);
00095 
00097   NFFT(vrand_shifted_unit_double)(p.x, p.M_total);
00098 
00100   if (p.flags & PRE_ONE_PSI)
00101     NFFT(precompute_one_psi)(&p);
00102 
00104   SOLVER(init_complex)(&ip, (NFFT(mv_plan_complex)*) (&p));
00105 
00107   NFFT(vrand_unit_complex)(ip.y, p.M_total);
00108   NFFT(vpr_complex)(ip.y, p.M_total, "Given data, vector y");
00109 
00111   for (k = 0; k < p.N_total; k++)
00112     ip.f_hat_iter[k] = NFFT_K(0.0);
00113 
00114   NFFT(vpr_complex)(ip.f_hat_iter, p.N_total,
00115       "Initial guess, vector f_hat_iter");
00116 
00118   error_str = NFFT(check)(&p);
00119   if (error_str != 0)
00120   {
00121     printf("Error in nfft module: %s\n", error_str);
00122     return;
00123   }
00124 
00125   NFFT_CSWAP(ip.f_hat_iter, p.f_hat);
00126   NFFT(trafo)(&p);
00127   NFFT(vpr_complex)(p.f, p.M_total, "Data fit, vector f");
00128   NFFT_CSWAP(ip.f_hat_iter, p.f_hat);
00129 
00130   SOLVER(before_loop_complex)(&ip);
00131   printf("\n Residual r=%" NFFT__FES__ "\n", ip.dot_r_iter);
00132 
00133   for (l = 0; l < iter; l++)
00134   {
00135     printf("\n********** Iteration l=%d **********\n", l);
00136     SOLVER(loop_one_step_complex)(&ip);
00137     NFFT(vpr_complex)(ip.f_hat_iter, p.N_total,
00138         "Approximate solution, vector f_hat_iter");
00139 
00140     NFFT_CSWAP(ip.f_hat_iter, p.f_hat);
00141     NFFT(trafo)(&p);
00142     NFFT(vpr_complex)(p.f, p.M_total, "Data fit, vector f");
00143     NFFT_CSWAP(ip.f_hat_iter, p.f_hat);
00144 
00145     printf("\n Residual r=%"  NFFT__FES__ "\n", ip.dot_r_iter);
00146   }
00147 
00148   SOLVER(finalize_complex)(&ip);
00149   NFFT(finalize)(&p);
00150 }
00151 
00153 int main(void)
00154 {
00155   printf("\n Computing a one dimensional inverse nfft\n");
00156 
00157   simple_test_solver_nfft_1d(16, 8, 9);
00158 
00159   return EXIT_SUCCESS;
00160 }
00161 /* \} */