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