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