![]() |
3.2.2 |
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 /* \} */