00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00026 #include <complex.h>
00027
00028 #include "nfft3util.h"
00029 #include "nfft3.h"
00030
00031 void solver_init_advanced_complex(solver_plan_complex* ths, mv_plan_complex *mv, unsigned flags)
00032 {
00033 ths->mv = mv;
00034 ths->flags = flags;
00035
00036 ths->y = (fftw_complex*)nfft_malloc(ths->mv->M_total*sizeof(fftw_complex));
00037 ths->r_iter = (fftw_complex*)nfft_malloc(ths->mv->M_total*sizeof(fftw_complex));
00038 ths->f_hat_iter = (fftw_complex*)nfft_malloc(ths->mv->N_total*sizeof(fftw_complex));
00039 ths->p_hat_iter = (fftw_complex*)nfft_malloc(ths->mv->N_total*sizeof(fftw_complex));
00040
00041 if(ths->flags & LANDWEBER)
00042 ths->z_hat_iter = ths->p_hat_iter;
00043
00044 if(ths->flags & STEEPEST_DESCENT)
00045 {
00046 ths->z_hat_iter = ths->p_hat_iter;
00047 ths->v_iter = (fftw_complex*)nfft_malloc(ths->mv->M_total*sizeof(fftw_complex));
00048 }
00049
00050 if(ths->flags & CGNR)
00051 {
00052 ths->z_hat_iter = (fftw_complex*)nfft_malloc(ths->mv->N_total*sizeof(fftw_complex));
00053 ths->v_iter = (fftw_complex*)nfft_malloc(ths->mv->M_total*sizeof(fftw_complex));
00054 }
00055
00056 if(ths->flags & CGNE)
00057 ths->z_hat_iter = ths->p_hat_iter;
00058
00059 if(ths->flags & PRECOMPUTE_WEIGHT)
00060 ths->w = (double*) nfft_malloc(ths->mv->M_total*sizeof(double));
00061
00062 if(ths->flags & PRECOMPUTE_DAMP)
00063 ths->w_hat = (double*) nfft_malloc(ths->mv->N_total*sizeof(double));
00064 }
00065
00066 void solver_init_complex(solver_plan_complex* ths, mv_plan_complex *mv)
00067 {
00068 solver_init_advanced_complex(ths, mv, CGNR);
00069 }
00070
00071 void solver_before_loop_complex(solver_plan_complex* ths)
00072 {
00073 nfft_cp_complex(ths->mv->f_hat, ths->f_hat_iter, ths->mv->N_total);
00074
00075 NFFT_SWAP_complex(ths->r_iter, ths->mv->f);
00076 ths->mv->mv_trafo(ths->mv);
00077 NFFT_SWAP_complex(ths->r_iter, ths->mv->f);
00078
00079 nfft_upd_axpy_complex(ths->r_iter, -1.0, ths->y, ths->mv->M_total);
00080
00081 if((!(ths->flags & LANDWEBER)) || (ths->flags & NORMS_FOR_LANDWEBER))
00082 {
00083 if(ths->flags & PRECOMPUTE_WEIGHT)
00084 ths->dot_r_iter = nfft_dot_w_complex(ths->r_iter, ths->w,
00085 ths->mv->M_total);
00086 else
00087 ths->dot_r_iter = nfft_dot_complex(ths->r_iter, ths->mv->M_total);
00088 }
00089
00090
00091 if(ths->flags & PRECOMPUTE_WEIGHT)
00092 nfft_cp_w_complex(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
00093 else
00094 nfft_cp_complex(ths->mv->f, ths->r_iter, ths->mv->M_total);
00095
00096 NFFT_SWAP_complex(ths->z_hat_iter, ths->mv->f_hat);
00097 ths->mv->mv_adjoint(ths->mv);
00098 NFFT_SWAP_complex(ths->z_hat_iter, ths->mv->f_hat);
00099
00100 if((!(ths->flags & LANDWEBER)) || (ths->flags & NORMS_FOR_LANDWEBER))
00101 {
00102 if(ths->flags & PRECOMPUTE_DAMP)
00103 ths->dot_z_hat_iter = nfft_dot_w_complex(ths->z_hat_iter, ths->w_hat,
00104 ths->mv->N_total);
00105 else
00106 ths->dot_z_hat_iter = nfft_dot_complex(ths->z_hat_iter,
00107 ths->mv->N_total);
00108 }
00109
00110 if(ths->flags & CGNE)
00111 ths->dot_p_hat_iter = ths->dot_z_hat_iter;
00112
00113 if(ths->flags & CGNR)
00114 nfft_cp_complex(ths->p_hat_iter, ths->z_hat_iter, ths->mv->N_total);
00115 }
00116
00118 void solver_loop_one_step_landweber_complex(solver_plan_complex* ths)
00119 {
00120 if(ths->flags & PRECOMPUTE_DAMP)
00121 nfft_upd_xpawy_complex(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
00122 ths->z_hat_iter, ths->mv->N_total);
00123 else
00124 nfft_upd_xpay_complex(ths->f_hat_iter, ths->alpha_iter, ths->z_hat_iter,
00125 ths->mv->N_total);
00126
00127
00128 nfft_cp_complex(ths->mv->f_hat, ths->f_hat_iter, ths->mv->N_total);
00129
00130 NFFT_SWAP_complex(ths->r_iter,ths->mv->f);
00131 ths->mv->mv_trafo(ths->mv);
00132 NFFT_SWAP_complex(ths->r_iter,ths->mv->f);
00133
00134 nfft_upd_axpy_complex(ths->r_iter, -1.0, ths->y, ths->mv->M_total);
00135
00136 if(ths->flags & NORMS_FOR_LANDWEBER)
00137 {
00138 if(ths->flags & PRECOMPUTE_WEIGHT)
00139 ths->dot_r_iter = nfft_dot_w_complex(ths->r_iter,ths->w,
00140 ths->mv->M_total);
00141 else
00142 ths->dot_r_iter = nfft_dot_complex(ths->r_iter, ths->mv->M_total);
00143 }
00144
00145
00146 if(ths->flags & PRECOMPUTE_WEIGHT)
00147 nfft_cp_w_complex(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
00148 else
00149 nfft_cp_complex(ths->mv->f, ths->r_iter, ths->mv->M_total);
00150
00151 NFFT_SWAP_complex(ths->z_hat_iter,ths->mv->f_hat);
00152 ths->mv->mv_adjoint(ths->mv);
00153 NFFT_SWAP_complex(ths->z_hat_iter,ths->mv->f_hat);
00154
00155 if(ths->flags & NORMS_FOR_LANDWEBER)
00156 {
00157 if(ths->flags & PRECOMPUTE_DAMP)
00158 ths->dot_z_hat_iter = nfft_dot_w_complex(ths->z_hat_iter, ths->w_hat,
00159 ths->mv->N_total);
00160 else
00161 ths->dot_z_hat_iter = nfft_dot_complex(ths->z_hat_iter,
00162 ths->mv->N_total);
00163 }
00164 }
00165
00167 void solver_loop_one_step_steepest_descent_complex(solver_plan_complex *ths)
00168 {
00169 if(ths->flags & PRECOMPUTE_DAMP)
00170 nfft_cp_w_complex(ths->mv->f_hat, ths->w_hat, ths->z_hat_iter,
00171 ths->mv->N_total);
00172 else
00173 nfft_cp_complex(ths->mv->f_hat, ths->z_hat_iter, ths->mv->N_total);
00174
00175 NFFT_SWAP_complex(ths->v_iter,ths->mv->f);
00176 ths->mv->mv_trafo(ths->mv);
00177 NFFT_SWAP_complex(ths->v_iter,ths->mv->f);
00178
00179 if(ths->flags & PRECOMPUTE_WEIGHT)
00180 ths->dot_v_iter = nfft_dot_w_complex(ths->v_iter,ths->w,ths->mv->M_total);
00181 else
00182 ths->dot_v_iter = nfft_dot_complex(ths->v_iter, ths->mv->M_total);
00183
00184
00185 ths->alpha_iter = ths->dot_z_hat_iter / ths->dot_v_iter;
00186
00187
00188 if(ths->flags & PRECOMPUTE_DAMP)
00189 nfft_upd_xpawy_complex(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
00190 ths->z_hat_iter, ths->mv->N_total);
00191 else
00192 nfft_upd_xpay_complex(ths->f_hat_iter, ths->alpha_iter, ths->z_hat_iter,
00193 ths->mv->N_total);
00194
00195
00196 nfft_upd_xpay_complex(ths->r_iter, -ths->alpha_iter, ths->v_iter,
00197 ths->mv->M_total);
00198
00199 if(ths->flags & PRECOMPUTE_WEIGHT)
00200 ths->dot_r_iter = nfft_dot_w_complex(ths->r_iter,ths->w,ths->mv->M_total);
00201 else
00202 ths->dot_r_iter = nfft_dot_complex(ths->r_iter, ths->mv->M_total);
00203
00204
00205 if(ths->flags & PRECOMPUTE_WEIGHT)
00206 nfft_cp_w_complex(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
00207 else
00208 nfft_cp_complex(ths->mv->f, ths->r_iter, ths->mv->M_total);
00209
00210 NFFT_SWAP_complex(ths->z_hat_iter,ths->mv->f_hat);
00211 ths->mv->mv_adjoint(ths->mv);
00212 NFFT_SWAP_complex(ths->z_hat_iter,ths->mv->f_hat);
00213
00214 if(ths->flags & PRECOMPUTE_DAMP)
00215 ths->dot_z_hat_iter = nfft_dot_w_complex(ths->z_hat_iter, ths->w_hat,
00216 ths->mv->N_total);
00217 else
00218 ths->dot_z_hat_iter = nfft_dot_complex(ths->z_hat_iter, ths->mv->N_total);
00219 }
00220
00222 void solver_loop_one_step_cgnr_complex(solver_plan_complex *ths)
00223 {
00224 if(ths->flags & PRECOMPUTE_DAMP)
00225 nfft_cp_w_complex(ths->mv->f_hat, ths->w_hat, ths->p_hat_iter,
00226 ths->mv->N_total);
00227 else
00228 nfft_cp_complex(ths->mv->f_hat, ths->p_hat_iter, ths->mv->N_total);
00229
00230 NFFT_SWAP_complex(ths->v_iter,ths->mv->f);
00231 ths->mv->mv_trafo(ths->mv);
00232 NFFT_SWAP_complex(ths->v_iter,ths->mv->f);
00233
00234 if(ths->flags & PRECOMPUTE_WEIGHT)
00235 ths->dot_v_iter = nfft_dot_w_complex(ths->v_iter,ths->w,ths->mv->M_total);
00236 else
00237 ths->dot_v_iter = nfft_dot_complex(ths->v_iter, ths->mv->M_total);
00238
00239
00240 ths->alpha_iter = ths->dot_z_hat_iter / ths->dot_v_iter;
00241
00242
00243 if(ths->flags & PRECOMPUTE_DAMP)
00244 nfft_upd_xpawy_complex(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
00245 ths->p_hat_iter, ths->mv->N_total);
00246 else
00247 nfft_upd_xpay_complex(ths->f_hat_iter, ths->alpha_iter, ths->p_hat_iter,
00248 ths->mv->N_total);
00249
00250
00251 nfft_upd_xpay_complex(ths->r_iter, -ths->alpha_iter, ths->v_iter,
00252 ths->mv->M_total);
00253
00254 if(ths->flags & PRECOMPUTE_WEIGHT)
00255 ths->dot_r_iter = nfft_dot_w_complex(ths->r_iter,ths->w,ths->mv->M_total);
00256 else
00257 ths->dot_r_iter = nfft_dot_complex(ths->r_iter, ths->mv->M_total);
00258
00259
00260 if(ths->flags & PRECOMPUTE_WEIGHT)
00261 nfft_cp_w_complex(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
00262 else
00263 nfft_cp_complex(ths->mv->f, ths->r_iter, ths->mv->M_total);
00264
00265 NFFT_SWAP_complex(ths->z_hat_iter,ths->mv->f_hat);
00266 ths->mv->mv_adjoint(ths->mv);
00267 NFFT_SWAP_complex(ths->z_hat_iter,ths->mv->f_hat);
00268
00269 ths->dot_z_hat_iter_old = ths->dot_z_hat_iter;
00270 if(ths->flags & PRECOMPUTE_DAMP)
00271 ths->dot_z_hat_iter = nfft_dot_w_complex(ths->z_hat_iter, ths->w_hat,
00272 ths->mv->N_total);
00273 else
00274 ths->dot_z_hat_iter = nfft_dot_complex(ths->z_hat_iter, ths->mv->N_total);
00275
00276
00277 ths->beta_iter = ths->dot_z_hat_iter / ths->dot_z_hat_iter_old;
00278
00279
00280 nfft_upd_axpy_complex(ths->p_hat_iter, ths->beta_iter, ths->z_hat_iter,
00281 ths->mv->N_total);
00282 }
00283
00285 void solver_loop_one_step_cgne_complex(solver_plan_complex *ths)
00286 {
00287 ths->alpha_iter = ths->dot_r_iter / ths->dot_p_hat_iter;
00288
00289
00290 if(ths->flags & PRECOMPUTE_DAMP)
00291 nfft_upd_xpawy_complex(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
00292 ths->p_hat_iter, ths->mv->N_total);
00293 else
00294 nfft_upd_xpay_complex(ths->f_hat_iter, ths->alpha_iter, ths->p_hat_iter,
00295 ths->mv->N_total);
00296
00297
00298 if(ths->flags & PRECOMPUTE_DAMP)
00299 nfft_cp_w_complex(ths->mv->f_hat, ths->w_hat, ths->p_hat_iter,
00300 ths->mv->N_total);
00301 else
00302 nfft_cp_complex(ths->mv->f_hat, ths->p_hat_iter, ths->mv->N_total);
00303
00304 ths->mv->mv_trafo(ths->mv);
00305
00306 nfft_upd_xpay_complex(ths->r_iter, -ths->alpha_iter, ths->mv->f,
00307 ths->mv->M_total);
00308
00309 ths->dot_r_iter_old = ths->dot_r_iter;
00310 if(ths->flags & PRECOMPUTE_WEIGHT)
00311 ths->dot_r_iter = nfft_dot_w_complex(ths->r_iter,ths->w,ths->mv->M_total);
00312 else
00313 ths->dot_r_iter = nfft_dot_complex(ths->r_iter, ths->mv->M_total);
00314
00315
00316 ths->beta_iter = ths->dot_r_iter / ths->dot_r_iter_old;
00317
00318
00319 if(ths->flags & PRECOMPUTE_WEIGHT)
00320 nfft_cp_w_complex(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
00321 else
00322 nfft_cp_complex(ths->mv->f, ths->r_iter, ths->mv->M_total);
00323
00324 ths->mv->mv_adjoint(ths->mv);
00325
00326 nfft_upd_axpy_complex(ths->p_hat_iter, ths->beta_iter, ths->mv->f_hat,
00327 ths->mv->N_total);
00328
00329 if(ths->flags & PRECOMPUTE_DAMP)
00330 ths->dot_p_hat_iter = nfft_dot_w_complex(ths->p_hat_iter, ths->w_hat,
00331 ths->mv->N_total);
00332 else
00333 ths->dot_p_hat_iter = nfft_dot_complex(ths->p_hat_iter, ths->mv->N_total);
00334 }
00335
00337 void solver_loop_one_step_complex(solver_plan_complex *ths)
00338 {
00339 if(ths->flags & LANDWEBER)
00340 solver_loop_one_step_landweber_complex(ths);
00341
00342 if(ths->flags & STEEPEST_DESCENT)
00343 solver_loop_one_step_steepest_descent_complex(ths);
00344
00345 if(ths->flags & CGNR)
00346 solver_loop_one_step_cgnr_complex(ths);
00347
00348 if(ths->flags & CGNE)
00349 solver_loop_one_step_cgne_complex(ths);
00350 }
00351
00353 void solver_finalize_complex(solver_plan_complex *ths)
00354 {
00355 if(ths->flags & PRECOMPUTE_WEIGHT)
00356 nfft_free(ths->w);
00357
00358 if(ths->flags & PRECOMPUTE_DAMP)
00359 nfft_free(ths->w_hat);
00360
00361 if(ths->flags & CGNR)
00362 {
00363 nfft_free(ths->v_iter);
00364 nfft_free(ths->z_hat_iter);
00365 }
00366
00367 if(ths->flags & STEEPEST_DESCENT)
00368 nfft_free(ths->v_iter);
00369
00370 nfft_free(ths->p_hat_iter);
00371 nfft_free(ths->f_hat_iter);
00372
00373 nfft_free(ths->r_iter);
00374 nfft_free(ths->y);
00375 }
00378
00379
00380
00381
00382 void solver_init_advanced_double(solver_plan_double* ths, mv_plan_double *mv, unsigned flags)
00383 {
00384 ths->mv = mv;
00385 ths->flags = flags;
00386
00387 ths->y = (double*)nfft_malloc(ths->mv->M_total*sizeof(double));
00388 ths->r_iter = (double*)nfft_malloc(ths->mv->M_total*sizeof(double));
00389 ths->f_hat_iter = (double*)nfft_malloc(ths->mv->N_total*sizeof(double));
00390 ths->p_hat_iter = (double*)nfft_malloc(ths->mv->N_total*sizeof(double));
00391
00392 if(ths->flags & LANDWEBER)
00393 ths->z_hat_iter = ths->p_hat_iter;
00394
00395 if(ths->flags & STEEPEST_DESCENT)
00396 {
00397 ths->z_hat_iter = ths->p_hat_iter;
00398 ths->v_iter = (double*)nfft_malloc(ths->mv->M_total*sizeof(double));
00399 }
00400
00401 if(ths->flags & CGNR)
00402 {
00403 ths->z_hat_iter = (double*)nfft_malloc(ths->mv->N_total*sizeof(double));
00404 ths->v_iter = (double*)nfft_malloc(ths->mv->M_total*sizeof(double));
00405 }
00406
00407 if(ths->flags & CGNE)
00408 ths->z_hat_iter = ths->p_hat_iter;
00409
00410 if(ths->flags & PRECOMPUTE_WEIGHT)
00411 ths->w = (double*) nfft_malloc(ths->mv->M_total*sizeof(double));
00412
00413 if(ths->flags & PRECOMPUTE_DAMP)
00414 ths->w_hat = (double*) nfft_malloc(ths->mv->N_total*sizeof(double));
00415 }
00416
00417 void solver_init_double(solver_plan_double* ths, mv_plan_double *mv)
00418 {
00419 solver_init_advanced_double(ths, mv, CGNR);
00420 }
00421
00422 void solver_before_loop_double(solver_plan_double* ths)
00423 {
00424 nfft_cp_double(ths->mv->f_hat, ths->f_hat_iter, ths->mv->N_total);
00425
00426 NFFT_SWAP_double(ths->r_iter, ths->mv->f);
00427 ths->mv->mv_trafo(ths->mv);
00428 NFFT_SWAP_double(ths->r_iter, ths->mv->f);
00429
00430 nfft_upd_axpy_double(ths->r_iter, -1.0, ths->y, ths->mv->M_total);
00431
00432 if((!(ths->flags & LANDWEBER)) || (ths->flags & NORMS_FOR_LANDWEBER))
00433 {
00434 if(ths->flags & PRECOMPUTE_WEIGHT)
00435 ths->dot_r_iter = nfft_dot_w_double(ths->r_iter, ths->w,
00436 ths->mv->M_total);
00437 else
00438 ths->dot_r_iter = nfft_dot_double(ths->r_iter, ths->mv->M_total);
00439 }
00440
00441
00442 if(ths->flags & PRECOMPUTE_WEIGHT)
00443 nfft_cp_w_double(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
00444 else
00445 nfft_cp_double(ths->mv->f, ths->r_iter, ths->mv->M_total);
00446
00447 NFFT_SWAP_double(ths->z_hat_iter, ths->mv->f_hat);
00448 ths->mv->mv_adjoint(ths->mv);
00449 NFFT_SWAP_double(ths->z_hat_iter, ths->mv->f_hat);
00450
00451 if((!(ths->flags & LANDWEBER)) || (ths->flags & NORMS_FOR_LANDWEBER))
00452 {
00453 if(ths->flags & PRECOMPUTE_DAMP)
00454 ths->dot_z_hat_iter = nfft_dot_w_double(ths->z_hat_iter, ths->w_hat,
00455 ths->mv->N_total);
00456 else
00457 ths->dot_z_hat_iter = nfft_dot_double(ths->z_hat_iter,
00458 ths->mv->N_total);
00459 }
00460
00461 if(ths->flags & CGNE)
00462 ths->dot_p_hat_iter = ths->dot_z_hat_iter;
00463
00464 if(ths->flags & CGNR)
00465 nfft_cp_double(ths->p_hat_iter, ths->z_hat_iter, ths->mv->N_total);
00466 }
00467
00469 void solver_loop_one_step_landweber_double(solver_plan_double* ths)
00470 {
00471 if(ths->flags & PRECOMPUTE_DAMP)
00472 nfft_upd_xpawy_double(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
00473 ths->z_hat_iter, ths->mv->N_total);
00474 else
00475 nfft_upd_xpay_double(ths->f_hat_iter, ths->alpha_iter, ths->z_hat_iter,
00476 ths->mv->N_total);
00477
00478
00479 nfft_cp_double(ths->mv->f_hat, ths->f_hat_iter, ths->mv->N_total);
00480
00481 NFFT_SWAP_double(ths->r_iter,ths->mv->f);
00482 ths->mv->mv_trafo(ths->mv);
00483 NFFT_SWAP_double(ths->r_iter,ths->mv->f);
00484
00485 nfft_upd_axpy_double(ths->r_iter, -1.0, ths->y, ths->mv->M_total);
00486
00487 if(ths->flags & NORMS_FOR_LANDWEBER)
00488 {
00489 if(ths->flags & PRECOMPUTE_WEIGHT)
00490 ths->dot_r_iter = nfft_dot_w_double(ths->r_iter,ths->w,
00491 ths->mv->M_total);
00492 else
00493 ths->dot_r_iter = nfft_dot_double(ths->r_iter, ths->mv->M_total);
00494 }
00495
00496
00497 if(ths->flags & PRECOMPUTE_WEIGHT)
00498 nfft_cp_w_double(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
00499 else
00500 nfft_cp_double(ths->mv->f, ths->r_iter, ths->mv->M_total);
00501
00502 NFFT_SWAP_double(ths->z_hat_iter,ths->mv->f_hat);
00503 ths->mv->mv_adjoint(ths->mv);
00504 NFFT_SWAP_double(ths->z_hat_iter,ths->mv->f_hat);
00505
00506 if(ths->flags & NORMS_FOR_LANDWEBER)
00507 {
00508 if(ths->flags & PRECOMPUTE_DAMP)
00509 ths->dot_z_hat_iter = nfft_dot_w_double(ths->z_hat_iter, ths->w_hat,
00510 ths->mv->N_total);
00511 else
00512 ths->dot_z_hat_iter = nfft_dot_double(ths->z_hat_iter,
00513 ths->mv->N_total);
00514 }
00515 }
00516
00518 void solver_loop_one_step_steepest_descent_double(solver_plan_double *ths)
00519 {
00520 if(ths->flags & PRECOMPUTE_DAMP)
00521 nfft_cp_w_double(ths->mv->f_hat, ths->w_hat, ths->z_hat_iter,
00522 ths->mv->N_total);
00523 else
00524 nfft_cp_double(ths->mv->f_hat, ths->z_hat_iter, ths->mv->N_total);
00525
00526 NFFT_SWAP_double(ths->v_iter,ths->mv->f);
00527 ths->mv->mv_trafo(ths->mv);
00528 NFFT_SWAP_double(ths->v_iter,ths->mv->f);
00529
00530 if(ths->flags & PRECOMPUTE_WEIGHT)
00531 ths->dot_v_iter = nfft_dot_w_double(ths->v_iter,ths->w,ths->mv->M_total);
00532 else
00533 ths->dot_v_iter = nfft_dot_double(ths->v_iter, ths->mv->M_total);
00534
00535
00536 ths->alpha_iter = ths->dot_z_hat_iter / ths->dot_v_iter;
00537
00538
00539 if(ths->flags & PRECOMPUTE_DAMP)
00540 nfft_upd_xpawy_double(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
00541 ths->z_hat_iter, ths->mv->N_total);
00542 else
00543 nfft_upd_xpay_double(ths->f_hat_iter, ths->alpha_iter, ths->z_hat_iter,
00544 ths->mv->N_total);
00545
00546
00547 nfft_upd_xpay_double(ths->r_iter, -ths->alpha_iter, ths->v_iter,
00548 ths->mv->M_total);
00549
00550 if(ths->flags & PRECOMPUTE_WEIGHT)
00551 ths->dot_r_iter = nfft_dot_w_double(ths->r_iter,ths->w,ths->mv->M_total);
00552 else
00553 ths->dot_r_iter = nfft_dot_double(ths->r_iter, ths->mv->M_total);
00554
00555
00556 if(ths->flags & PRECOMPUTE_WEIGHT)
00557 nfft_cp_w_double(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
00558 else
00559 nfft_cp_double(ths->mv->f, ths->r_iter, ths->mv->M_total);
00560
00561 NFFT_SWAP_double(ths->z_hat_iter,ths->mv->f_hat);
00562 ths->mv->mv_adjoint(ths->mv);
00563 NFFT_SWAP_double(ths->z_hat_iter,ths->mv->f_hat);
00564
00565 if(ths->flags & PRECOMPUTE_DAMP)
00566 ths->dot_z_hat_iter = nfft_dot_w_double(ths->z_hat_iter, ths->w_hat,
00567 ths->mv->N_total);
00568 else
00569 ths->dot_z_hat_iter = nfft_dot_double(ths->z_hat_iter, ths->mv->N_total);
00570 }
00571
00573 void solver_loop_one_step_cgnr_double(solver_plan_double *ths)
00574 {
00575 if(ths->flags & PRECOMPUTE_DAMP)
00576 nfft_cp_w_double(ths->mv->f_hat, ths->w_hat, ths->p_hat_iter,
00577 ths->mv->N_total);
00578 else
00579 nfft_cp_double(ths->mv->f_hat, ths->p_hat_iter, ths->mv->N_total);
00580
00581 NFFT_SWAP_double(ths->v_iter,ths->mv->f);
00582 ths->mv->mv_trafo(ths->mv);
00583 NFFT_SWAP_double(ths->v_iter,ths->mv->f);
00584
00585 if(ths->flags & PRECOMPUTE_WEIGHT)
00586 ths->dot_v_iter = nfft_dot_w_double(ths->v_iter,ths->w,ths->mv->M_total);
00587 else
00588 ths->dot_v_iter = nfft_dot_double(ths->v_iter, ths->mv->M_total);
00589
00590
00591 ths->alpha_iter = ths->dot_z_hat_iter / ths->dot_v_iter;
00592
00593
00594 if(ths->flags & PRECOMPUTE_DAMP)
00595 nfft_upd_xpawy_double(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
00596 ths->p_hat_iter, ths->mv->N_total);
00597 else
00598 nfft_upd_xpay_double(ths->f_hat_iter, ths->alpha_iter, ths->p_hat_iter,
00599 ths->mv->N_total);
00600
00601
00602 nfft_upd_xpay_double(ths->r_iter, -ths->alpha_iter, ths->v_iter,
00603 ths->mv->M_total);
00604
00605 if(ths->flags & PRECOMPUTE_WEIGHT)
00606 ths->dot_r_iter = nfft_dot_w_double(ths->r_iter,ths->w,ths->mv->M_total);
00607 else
00608 ths->dot_r_iter = nfft_dot_double(ths->r_iter, ths->mv->M_total);
00609
00610
00611 if(ths->flags & PRECOMPUTE_WEIGHT)
00612 nfft_cp_w_double(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
00613 else
00614 nfft_cp_double(ths->mv->f, ths->r_iter, ths->mv->M_total);
00615
00616 NFFT_SWAP_double(ths->z_hat_iter,ths->mv->f_hat);
00617 ths->mv->mv_adjoint(ths->mv);
00618 NFFT_SWAP_double(ths->z_hat_iter,ths->mv->f_hat);
00619
00620 ths->dot_z_hat_iter_old = ths->dot_z_hat_iter;
00621 if(ths->flags & PRECOMPUTE_DAMP)
00622 ths->dot_z_hat_iter = nfft_dot_w_double(ths->z_hat_iter, ths->w_hat,
00623 ths->mv->N_total);
00624 else
00625 ths->dot_z_hat_iter = nfft_dot_double(ths->z_hat_iter, ths->mv->N_total);
00626
00627
00628 ths->beta_iter = ths->dot_z_hat_iter / ths->dot_z_hat_iter_old;
00629
00630
00631 nfft_upd_axpy_double(ths->p_hat_iter, ths->beta_iter, ths->z_hat_iter,
00632 ths->mv->N_total);
00633 }
00634
00636 void solver_loop_one_step_cgne_double(solver_plan_double *ths)
00637 {
00638 ths->alpha_iter = ths->dot_r_iter / ths->dot_p_hat_iter;
00639
00640
00641 if(ths->flags & PRECOMPUTE_DAMP)
00642 nfft_upd_xpawy_double(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
00643 ths->p_hat_iter, ths->mv->N_total);
00644 else
00645 nfft_upd_xpay_double(ths->f_hat_iter, ths->alpha_iter, ths->p_hat_iter,
00646 ths->mv->N_total);
00647
00648
00649 if(ths->flags & PRECOMPUTE_DAMP)
00650 nfft_cp_w_double(ths->mv->f_hat, ths->w_hat, ths->p_hat_iter,
00651 ths->mv->N_total);
00652 else
00653 nfft_cp_double(ths->mv->f_hat, ths->p_hat_iter, ths->mv->N_total);
00654
00655 ths->mv->mv_trafo(ths->mv);
00656
00657 nfft_upd_xpay_double(ths->r_iter, -ths->alpha_iter, ths->mv->f,
00658 ths->mv->M_total);
00659
00660 ths->dot_r_iter_old = ths->dot_r_iter;
00661 if(ths->flags & PRECOMPUTE_WEIGHT)
00662 ths->dot_r_iter = nfft_dot_w_double(ths->r_iter,ths->w,ths->mv->M_total);
00663 else
00664 ths->dot_r_iter = nfft_dot_double(ths->r_iter, ths->mv->M_total);
00665
00666
00667 ths->beta_iter = ths->dot_r_iter / ths->dot_r_iter_old;
00668
00669
00670 if(ths->flags & PRECOMPUTE_WEIGHT)
00671 nfft_cp_w_double(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
00672 else
00673 nfft_cp_double(ths->mv->f, ths->r_iter, ths->mv->M_total);
00674
00675 ths->mv->mv_adjoint(ths->mv);
00676
00677 nfft_upd_axpy_double(ths->p_hat_iter, ths->beta_iter, ths->mv->f_hat,
00678 ths->mv->N_total);
00679
00680 if(ths->flags & PRECOMPUTE_DAMP)
00681 ths->dot_p_hat_iter = nfft_dot_w_double(ths->p_hat_iter, ths->w_hat,
00682 ths->mv->N_total);
00683 else
00684 ths->dot_p_hat_iter = nfft_dot_double(ths->p_hat_iter, ths->mv->N_total);
00685 }
00686
00688 void solver_loop_one_step_double(solver_plan_double *ths)
00689 {
00690 if(ths->flags & LANDWEBER)
00691 solver_loop_one_step_landweber_double(ths);
00692
00693 if(ths->flags & STEEPEST_DESCENT)
00694 solver_loop_one_step_steepest_descent_double(ths);
00695
00696 if(ths->flags & CGNR)
00697 solver_loop_one_step_cgnr_double(ths);
00698
00699 if(ths->flags & CGNE)
00700 solver_loop_one_step_cgne_double(ths);
00701 }
00702
00704 void solver_finalize_double(solver_plan_double *ths)
00705 {
00706 if(ths->flags & PRECOMPUTE_WEIGHT)
00707 nfft_free(ths->w);
00708
00709 if(ths->flags & PRECOMPUTE_DAMP)
00710 nfft_free(ths->w_hat);
00711
00712 if(ths->flags & CGNR)
00713 {
00714 nfft_free(ths->v_iter);
00715 nfft_free(ths->z_hat_iter);
00716 }
00717
00718 if(ths->flags & STEEPEST_DESCENT)
00719 nfft_free(ths->v_iter);
00720
00721 nfft_free(ths->p_hat_iter);
00722 nfft_free(ths->f_hat_iter);
00723
00724 nfft_free(ths->r_iter);
00725 nfft_free(ths->y);
00726 }