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