NFFT Logo 3.2.2
solver.c
Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 2002, 2012 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: solver.c 3896 2012-10-10 12:19:26Z tovo $ */
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 } /* void solver_before_loop */
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 } /* void solver_loop_one_step_landweber */
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 } /* void solver_loop_one_step_steepest_descent */
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 } /* void solver_loop_one_step_cgnr */
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 } /* void solver_loop_one_step_cgne */
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 } /* void solver_loop_one_step */
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 } /* void solver_before_loop */
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 } /* void solver_loop_one_step_landweber */
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 } /* void solver_loop_one_step_steepest_descent */
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 } /* void solver_loop_one_step_cgnr */
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 } /* void solver_loop_one_step_cgne */
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 } /* void solver_loop_one_step */
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 } 

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409