00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
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 #include "infft.h"
00033
00042 static double my_weight(double z,double a,double b,double c)
00043 {
00044 return pow(0.25-z*z,b)/(c+pow(fabs(z),2*a));
00045 }
00046
00048 static void glacier(int N,int M)
00049 {
00050 int j,k,k0,k1,l,my_N[2],my_n[2];
00051 double tmp_y;
00052 nfft_plan p;
00053 solver_plan_complex ip;
00054 FILE* fp;
00055
00056
00057 my_N[0]=N; my_n[0]=X(next_power_of_2)(N);
00058 my_N[1]=N; my_n[1]=X(next_power_of_2)(N);
00059 nfft_init_guru(&p, 2, my_N, M, my_n, 6,
00060 PRE_PHI_HUT| PRE_FULL_PSI|
00061 MALLOC_X| MALLOC_F_HAT| MALLOC_F|
00062 FFTW_INIT| FFT_OUT_OF_PLACE,
00063 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00064
00065
00066 solver_init_advanced_complex(&ip,(nfft_mv_plan_complex*)(&p), CGNE| PRECOMPUTE_DAMP);
00067 fprintf(stderr,"Using the generic solver!");
00068
00069
00070 fp=fopen("input_data.dat","r");
00071 for(j=0;j<p.M_total;j++)
00072 {
00073 fscanf(fp,"%le %le %le",&p.x[2*j+0],&p.x[2*j+1],&tmp_y);
00074 ip.y[j]=tmp_y;
00075 }
00076 fclose(fp);
00077
00078
00079 if(p.nfft_flags & PRE_ONE_PSI)
00080 nfft_precompute_one_psi(&p);
00081
00082
00083 if(ip.flags & PRECOMPUTE_DAMP)
00084 for(k0=0;k0<p.N[0];k0++)
00085 for(k1=0;k1<p.N[1];k1++)
00086 ip.w_hat[k0*p.N[1]+k1]=
00087 my_weight(((double)(k0-p.N[0]/2))/p.N[0],0.5,3,0.001)*
00088 my_weight(((double)(k1-p.N[1]/2))/p.N[1],0.5,3,0.001);
00089
00090
00091 for(k=0;k<p.N_total;k++)
00092 ip.f_hat_iter[k]=0;
00093
00094
00095 solver_before_loop_complex(&ip);
00096 for(l=0;l<40;l++)
00097 {
00098 fprintf(stderr,"Residual ||r||=%e,\n",sqrt(ip.dot_r_iter));
00099 solver_loop_one_step_complex(&ip);
00100 }
00101
00102 for(k=0;k<p.N_total;k++)
00103 printf("%le %le\n",creal(ip.f_hat_iter[k]),cimag(ip.f_hat_iter[k]));
00104
00105 solver_finalize_complex(&ip);
00106 nfft_finalize(&p);
00107 }
00108
00110 static void glacier_cv(int N,int M,int M_cv,unsigned solver_flags)
00111 {
00112 int j,k,k0,k1,l,my_N[2],my_n[2];
00113 double tmp_y,r;
00114 nfft_plan p,cp;
00115 solver_plan_complex ip;
00116 double _Complex* cp_y;
00117 FILE* fp;
00118 int M_re=M-M_cv;
00119
00120
00121 my_N[0]=N; my_n[0]=X(next_power_of_2)(N);
00122 my_N[1]=N; my_n[1]=X(next_power_of_2)(N);
00123 nfft_init_guru(&p, 2, my_N, M_re, my_n, 6,
00124 PRE_PHI_HUT| PRE_FULL_PSI|
00125 MALLOC_X| MALLOC_F_HAT| MALLOC_F|
00126 FFTW_INIT| FFT_OUT_OF_PLACE,
00127 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00128
00129
00130
00131 solver_init_advanced_complex(&ip,(nfft_mv_plan_complex*)(&p), solver_flags);
00132
00133
00134 cp_y = (double _Complex*) nfft_malloc(M*sizeof(double _Complex));
00135 nfft_init_guru(&cp, 2, my_N, M, my_n, 6,
00136 PRE_PHI_HUT| PRE_FULL_PSI|
00137 MALLOC_X| MALLOC_F|
00138 FFTW_INIT| FFT_OUT_OF_PLACE,
00139 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00140
00141 cp.f_hat=ip.f_hat_iter;
00142
00143
00144 fp=fopen("input_data.dat","r");
00145 for(j=0;j<cp.M_total;j++)
00146 {
00147 fscanf(fp,"%le %le %le",&cp.x[2*j+0],&cp.x[2*j+1],&tmp_y);
00148 cp_y[j]=tmp_y;
00149 }
00150 fclose(fp);
00151
00152
00153 for(j=0;j<p.M_total;j++)
00154 {
00155 p.x[2*j+0]=cp.x[2*j+0];
00156 p.x[2*j+1]=cp.x[2*j+1];
00157 ip.y[j]=tmp_y;
00158 }
00159
00160
00161 if(p.nfft_flags & PRE_ONE_PSI)
00162 nfft_precompute_one_psi(&p);
00163
00164
00165 if(cp.nfft_flags & PRE_ONE_PSI)
00166 nfft_precompute_one_psi(&cp);
00167
00168
00169 if(ip.flags & PRECOMPUTE_DAMP)
00170 for(k0=0;k0<p.N[0];k0++)
00171 for(k1=0;k1<p.N[1];k1++)
00172 ip.w_hat[k0*p.N[1]+k1]=
00173 my_weight(((double)(k0-p.N[0]/2))/p.N[0],0.5,3,0.001)*
00174 my_weight(((double)(k1-p.N[1]/2))/p.N[1],0.5,3,0.001);
00175
00176
00177 for(k=0;k<p.N_total;k++)
00178 ip.f_hat_iter[k]=0;
00179
00180
00181 solver_before_loop_complex(&ip);
00182
00183 for(l=0;l<40;l++)
00184 solver_loop_one_step_complex(&ip);
00185
00186
00187
00188 NFFT_SWAP_complex(p.f_hat,ip.f_hat_iter);
00189 nfft_trafo(&p);
00190 NFFT_SWAP_complex(p.f_hat,ip.f_hat_iter);
00191 nfft_upd_axpy_complex(p.f,-1,ip.y,M_re);
00192 r=sqrt(nfft_dot_complex(p.f,M_re)/nfft_dot_complex(cp_y,M));
00193 fprintf(stderr,"r=%1.2e, ",r);
00194 printf("$%1.1e$ & ",r);
00195
00196 nfft_trafo(&cp);
00197 nfft_upd_axpy_complex(&cp.f[M_re],-1,&cp_y[M_re],M_cv);
00198 r=sqrt(nfft_dot_complex(&cp.f[M_re],M_cv)/nfft_dot_complex(cp_y,M));
00199 fprintf(stderr,"r_1=%1.2e\t",r);
00200 printf("$%1.1e$ & ",r);
00201
00202 nfft_finalize(&cp);
00203 solver_finalize_complex(&ip);
00204 nfft_finalize(&p);
00205 }
00206
00207
00209 int main(int argc, char **argv)
00210 {
00211 int M_cv;
00212
00213 if(argc<3)
00214 {
00215 fprintf(stderr,"Call this program from the Matlab script glacier.m!");
00216 exit(-1);
00217 }
00218
00219 if(argc==3)
00220 glacier(atoi(argv[1]),atoi(argv[2]));
00221 else
00222 for(M_cv=atoi(argv[3]);M_cv<=atoi(argv[5]);M_cv+=atoi(argv[4]))
00223 {
00224 fprintf(stderr,"\nM_cv=%d,\t",M_cv);
00225 printf("$%d$ & ",M_cv);
00226 fprintf(stderr,"cgne+damp: ");
00227 glacier_cv(atoi(argv[1]),atoi(argv[2]),M_cv,CGNE| PRECOMPUTE_DAMP);
00228
00229
00230 fprintf(stderr,"cgnr: ");
00231 glacier_cv(atoi(argv[1]),atoi(argv[2]),M_cv,CGNR);
00232 fprintf(stderr,"cgnr: ");
00233 glacier_cv(atoi(argv[1])/4,atoi(argv[2]),M_cv,CGNR);
00234 printf("XXX \\\\\n");
00235 }
00236
00237 fprintf(stderr,"\n");
00238
00239 return 1;
00240 }
00241