NFFT Logo 3.2.3
solver.c
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2002, 2012 Jens Keiner, Stefan Kunis, Daniel Potts
3  *
4  * This program is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU General Public License as published by the Free Software
6  * Foundation; either version 2 of the License, or (at your option) any later
7  * version.
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU General Public License along with
15  * this program; if not, write to the Free Software Foundation, Inc., 51
16  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18 
19 /* $Id: solver.c 3896 2012-10-10 12:19:26Z tovo $ */
20 
25 #include "config.h"
26 
27 #ifdef HAVE_COMPLEX_H
28 #include <complex.h>
29 #endif
30 #include "nfft3util.h"
31 #include "nfft3.h"
32 
33 void solver_init_advanced_complex(solver_plan_complex* ths, nfft_mv_plan_complex *mv, unsigned flags)
34 {
35  ths->mv = mv;
36  ths->flags = flags;
37 
38  ths->y = (fftw_complex*)nfft_malloc(ths->mv->M_total*sizeof(fftw_complex));
39  ths->r_iter = (fftw_complex*)nfft_malloc(ths->mv->M_total*sizeof(fftw_complex));
40  ths->f_hat_iter = (fftw_complex*)nfft_malloc(ths->mv->N_total*sizeof(fftw_complex));
41  ths->p_hat_iter = (fftw_complex*)nfft_malloc(ths->mv->N_total*sizeof(fftw_complex));
42 
43  if(ths->flags & LANDWEBER)
44  ths->z_hat_iter = ths->p_hat_iter;
45 
46  if(ths->flags & STEEPEST_DESCENT)
47  {
48  ths->z_hat_iter = ths->p_hat_iter;
49  ths->v_iter = (fftw_complex*)nfft_malloc(ths->mv->M_total*sizeof(fftw_complex));
50  }
51 
52  if(ths->flags & CGNR)
53  {
54  ths->z_hat_iter = (fftw_complex*)nfft_malloc(ths->mv->N_total*sizeof(fftw_complex));
55  ths->v_iter = (fftw_complex*)nfft_malloc(ths->mv->M_total*sizeof(fftw_complex));
56  }
57 
58  if(ths->flags & CGNE)
59  ths->z_hat_iter = ths->p_hat_iter;
60 
61  if(ths->flags & PRECOMPUTE_WEIGHT)
62  ths->w = (double*) nfft_malloc(ths->mv->M_total*sizeof(double));
63 
64  if(ths->flags & PRECOMPUTE_DAMP)
65  ths->w_hat = (double*) nfft_malloc(ths->mv->N_total*sizeof(double));
66 }
67 
68 void solver_init_complex(solver_plan_complex* ths, nfft_mv_plan_complex *mv)
69 {
70  solver_init_advanced_complex(ths, mv, CGNR);
71 }
72 
73 void solver_before_loop_complex(solver_plan_complex* ths)
74 {
75  nfft_cp_complex(ths->mv->f_hat, ths->f_hat_iter, ths->mv->N_total);
76 
77  NFFT_SWAP_complex(ths->r_iter, ths->mv->f);
78  ths->mv->mv_trafo(ths->mv);
79  NFFT_SWAP_complex(ths->r_iter, ths->mv->f);
80 
81  nfft_upd_axpy_complex(ths->r_iter, -1.0, ths->y, ths->mv->M_total);
82 
83  if((!(ths->flags & LANDWEBER)) || (ths->flags & NORMS_FOR_LANDWEBER))
84  {
85  if(ths->flags & PRECOMPUTE_WEIGHT)
86  ths->dot_r_iter = nfft_dot_w_complex(ths->r_iter, ths->w,
87  ths->mv->M_total);
88  else
89  ths->dot_r_iter = nfft_dot_complex(ths->r_iter, ths->mv->M_total);
90  }
91 
92  /*-----------------*/
93  if(ths->flags & PRECOMPUTE_WEIGHT)
94  nfft_cp_w_complex(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
95  else
96  nfft_cp_complex(ths->mv->f, ths->r_iter, ths->mv->M_total);
97 
98  NFFT_SWAP_complex(ths->z_hat_iter, ths->mv->f_hat);
99  ths->mv->mv_adjoint(ths->mv);
100  NFFT_SWAP_complex(ths->z_hat_iter, ths->mv->f_hat);
101 
102  if((!(ths->flags & LANDWEBER)) || (ths->flags & NORMS_FOR_LANDWEBER))
103  {
104  if(ths->flags & PRECOMPUTE_DAMP)
105  ths->dot_z_hat_iter = nfft_dot_w_complex(ths->z_hat_iter, ths->w_hat,
106  ths->mv->N_total);
107  else
108  ths->dot_z_hat_iter = nfft_dot_complex(ths->z_hat_iter,
109  ths->mv->N_total);
110  }
111 
112  if(ths->flags & CGNE)
113  ths->dot_p_hat_iter = ths->dot_z_hat_iter;
114 
115  if(ths->flags & CGNR)
116  nfft_cp_complex(ths->p_hat_iter, ths->z_hat_iter, ths->mv->N_total);
117 } /* void solver_before_loop */
118 
120 static void solver_loop_one_step_landweber_complex(solver_plan_complex* ths)
121 {
122  if(ths->flags & PRECOMPUTE_DAMP)
123  nfft_upd_xpawy_complex(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
124  ths->z_hat_iter, ths->mv->N_total);
125  else
126  nfft_upd_xpay_complex(ths->f_hat_iter, ths->alpha_iter, ths->z_hat_iter,
127  ths->mv->N_total);
128 
129  /*-----------------*/
130  nfft_cp_complex(ths->mv->f_hat, ths->f_hat_iter, ths->mv->N_total);
131 
132  NFFT_SWAP_complex(ths->r_iter,ths->mv->f);
133  ths->mv->mv_trafo(ths->mv);
134  NFFT_SWAP_complex(ths->r_iter,ths->mv->f);
135 
136  nfft_upd_axpy_complex(ths->r_iter, -1.0, ths->y, ths->mv->M_total);
137 
138  if(ths->flags & NORMS_FOR_LANDWEBER)
139  {
140  if(ths->flags & PRECOMPUTE_WEIGHT)
141  ths->dot_r_iter = nfft_dot_w_complex(ths->r_iter,ths->w,
142  ths->mv->M_total);
143  else
144  ths->dot_r_iter = nfft_dot_complex(ths->r_iter, ths->mv->M_total);
145  }
146 
147  /*-----------------*/
148  if(ths->flags & PRECOMPUTE_WEIGHT)
149  nfft_cp_w_complex(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
150  else
151  nfft_cp_complex(ths->mv->f, ths->r_iter, ths->mv->M_total);
152 
153  NFFT_SWAP_complex(ths->z_hat_iter,ths->mv->f_hat);
154  ths->mv->mv_adjoint(ths->mv);
155  NFFT_SWAP_complex(ths->z_hat_iter,ths->mv->f_hat);
156 
157  if(ths->flags & NORMS_FOR_LANDWEBER)
158  {
159  if(ths->flags & PRECOMPUTE_DAMP)
160  ths->dot_z_hat_iter = nfft_dot_w_complex(ths->z_hat_iter, ths->w_hat,
161  ths->mv->N_total);
162  else
163  ths->dot_z_hat_iter = nfft_dot_complex(ths->z_hat_iter,
164  ths->mv->N_total);
165  }
166 } /* void solver_loop_one_step_landweber */
167 
169 static void solver_loop_one_step_steepest_descent_complex(solver_plan_complex *ths)
170 {
171  if(ths->flags & PRECOMPUTE_DAMP)
172  nfft_cp_w_complex(ths->mv->f_hat, ths->w_hat, ths->z_hat_iter,
173  ths->mv->N_total);
174  else
175  nfft_cp_complex(ths->mv->f_hat, ths->z_hat_iter, ths->mv->N_total);
176 
177  NFFT_SWAP_complex(ths->v_iter,ths->mv->f);
178  ths->mv->mv_trafo(ths->mv);
179  NFFT_SWAP_complex(ths->v_iter,ths->mv->f);
180 
181  if(ths->flags & PRECOMPUTE_WEIGHT)
182  ths->dot_v_iter = nfft_dot_w_complex(ths->v_iter,ths->w,ths->mv->M_total);
183  else
184  ths->dot_v_iter = nfft_dot_complex(ths->v_iter, ths->mv->M_total);
185 
186  /*-----------------*/
187  ths->alpha_iter = ths->dot_z_hat_iter / ths->dot_v_iter;
188 
189  /*-----------------*/
190  if(ths->flags & PRECOMPUTE_DAMP)
191  nfft_upd_xpawy_complex(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
192  ths->z_hat_iter, ths->mv->N_total);
193  else
194  nfft_upd_xpay_complex(ths->f_hat_iter, ths->alpha_iter, ths->z_hat_iter,
195  ths->mv->N_total);
196 
197  /*-----------------*/
198  nfft_upd_xpay_complex(ths->r_iter, -ths->alpha_iter, ths->v_iter,
199  ths->mv->M_total);
200 
201  if(ths->flags & PRECOMPUTE_WEIGHT)
202  ths->dot_r_iter = nfft_dot_w_complex(ths->r_iter,ths->w,ths->mv->M_total);
203  else
204  ths->dot_r_iter = nfft_dot_complex(ths->r_iter, ths->mv->M_total);
205 
206  /*-----------------*/
207  if(ths->flags & PRECOMPUTE_WEIGHT)
208  nfft_cp_w_complex(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
209  else
210  nfft_cp_complex(ths->mv->f, ths->r_iter, ths->mv->M_total);
211 
212  NFFT_SWAP_complex(ths->z_hat_iter,ths->mv->f_hat);
213  ths->mv->mv_adjoint(ths->mv);
214  NFFT_SWAP_complex(ths->z_hat_iter,ths->mv->f_hat);
215 
216  if(ths->flags & PRECOMPUTE_DAMP)
217  ths->dot_z_hat_iter = nfft_dot_w_complex(ths->z_hat_iter, ths->w_hat,
218  ths->mv->N_total);
219  else
220  ths->dot_z_hat_iter = nfft_dot_complex(ths->z_hat_iter, ths->mv->N_total);
221 } /* void solver_loop_one_step_steepest_descent */
222 
224 static void solver_loop_one_step_cgnr_complex(solver_plan_complex *ths)
225 {
226  if(ths->flags & PRECOMPUTE_DAMP)
227  nfft_cp_w_complex(ths->mv->f_hat, ths->w_hat, ths->p_hat_iter,
228  ths->mv->N_total);
229  else
230  nfft_cp_complex(ths->mv->f_hat, ths->p_hat_iter, ths->mv->N_total);
231 
232  NFFT_SWAP_complex(ths->v_iter,ths->mv->f);
233  ths->mv->mv_trafo(ths->mv);
234  NFFT_SWAP_complex(ths->v_iter,ths->mv->f);
235 
236  if(ths->flags & PRECOMPUTE_WEIGHT)
237  ths->dot_v_iter = nfft_dot_w_complex(ths->v_iter,ths->w,ths->mv->M_total);
238  else
239  ths->dot_v_iter = nfft_dot_complex(ths->v_iter, ths->mv->M_total);
240 
241  /*-----------------*/
242  ths->alpha_iter = ths->dot_z_hat_iter / ths->dot_v_iter;
243 
244  /*-----------------*/
245  if(ths->flags & PRECOMPUTE_DAMP)
246  nfft_upd_xpawy_complex(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
247  ths->p_hat_iter, ths->mv->N_total);
248  else
249  nfft_upd_xpay_complex(ths->f_hat_iter, ths->alpha_iter, ths->p_hat_iter,
250  ths->mv->N_total);
251 
252  /*-----------------*/
253  nfft_upd_xpay_complex(ths->r_iter, -ths->alpha_iter, ths->v_iter,
254  ths->mv->M_total);
255 
256  if(ths->flags & PRECOMPUTE_WEIGHT)
257  ths->dot_r_iter = nfft_dot_w_complex(ths->r_iter,ths->w,ths->mv->M_total);
258  else
259  ths->dot_r_iter = nfft_dot_complex(ths->r_iter, ths->mv->M_total);
260 
261  /*-----------------*/
262  if(ths->flags & PRECOMPUTE_WEIGHT)
263  nfft_cp_w_complex(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
264  else
265  nfft_cp_complex(ths->mv->f, ths->r_iter, ths->mv->M_total);
266 
267  NFFT_SWAP_complex(ths->z_hat_iter,ths->mv->f_hat);
268  ths->mv->mv_adjoint(ths->mv);
269  NFFT_SWAP_complex(ths->z_hat_iter,ths->mv->f_hat);
270 
271  ths->dot_z_hat_iter_old = ths->dot_z_hat_iter;
272  if(ths->flags & PRECOMPUTE_DAMP)
273  ths->dot_z_hat_iter = nfft_dot_w_complex(ths->z_hat_iter, ths->w_hat,
274  ths->mv->N_total);
275  else
276  ths->dot_z_hat_iter = nfft_dot_complex(ths->z_hat_iter, ths->mv->N_total);
277 
278  /*-----------------*/
279  ths->beta_iter = ths->dot_z_hat_iter / ths->dot_z_hat_iter_old;
280 
281  /*-----------------*/
282  nfft_upd_axpy_complex(ths->p_hat_iter, ths->beta_iter, ths->z_hat_iter,
283  ths->mv->N_total);
284 } /* void solver_loop_one_step_cgnr */
285 
287 static void solver_loop_one_step_cgne_complex(solver_plan_complex *ths)
288 {
289  ths->alpha_iter = ths->dot_r_iter / ths->dot_p_hat_iter;
290 
291  /*-----------------*/
292  if(ths->flags & PRECOMPUTE_DAMP)
293  nfft_upd_xpawy_complex(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
294  ths->p_hat_iter, ths->mv->N_total);
295  else
296  nfft_upd_xpay_complex(ths->f_hat_iter, ths->alpha_iter, ths->p_hat_iter,
297  ths->mv->N_total);
298 
299  /*-----------------*/
300  if(ths->flags & PRECOMPUTE_DAMP)
301  nfft_cp_w_complex(ths->mv->f_hat, ths->w_hat, ths->p_hat_iter,
302  ths->mv->N_total);
303  else
304  nfft_cp_complex(ths->mv->f_hat, ths->p_hat_iter, ths->mv->N_total);
305 
306  ths->mv->mv_trafo(ths->mv);
307 
308  nfft_upd_xpay_complex(ths->r_iter, -ths->alpha_iter, ths->mv->f,
309  ths->mv->M_total);
310 
311  ths->dot_r_iter_old = ths->dot_r_iter;
312  if(ths->flags & PRECOMPUTE_WEIGHT)
313  ths->dot_r_iter = nfft_dot_w_complex(ths->r_iter,ths->w,ths->mv->M_total);
314  else
315  ths->dot_r_iter = nfft_dot_complex(ths->r_iter, ths->mv->M_total);
316 
317  /*-----------------*/
318  ths->beta_iter = ths->dot_r_iter / ths->dot_r_iter_old;
319 
320  /*-----------------*/
321  if(ths->flags & PRECOMPUTE_WEIGHT)
322  nfft_cp_w_complex(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
323  else
324  nfft_cp_complex(ths->mv->f, ths->r_iter, ths->mv->M_total);
325 
326  ths->mv->mv_adjoint(ths->mv);
327 
328  nfft_upd_axpy_complex(ths->p_hat_iter, ths->beta_iter, ths->mv->f_hat,
329  ths->mv->N_total);
330 
331  if(ths->flags & PRECOMPUTE_DAMP)
332  ths->dot_p_hat_iter = nfft_dot_w_complex(ths->p_hat_iter, ths->w_hat,
333  ths->mv->N_total);
334  else
335  ths->dot_p_hat_iter = nfft_dot_complex(ths->p_hat_iter, ths->mv->N_total);
336 } /* void solver_loop_one_step_cgne */
337 
339 void solver_loop_one_step_complex(solver_plan_complex *ths)
340 {
341  if(ths->flags & LANDWEBER)
342  solver_loop_one_step_landweber_complex(ths);
343 
344  if(ths->flags & STEEPEST_DESCENT)
345  solver_loop_one_step_steepest_descent_complex(ths);
346 
347  if(ths->flags & CGNR)
348  solver_loop_one_step_cgnr_complex(ths);
349 
350  if(ths->flags & CGNE)
351  solver_loop_one_step_cgne_complex(ths);
352 } /* void solver_loop_one_step */
353 
355 void solver_finalize_complex(solver_plan_complex *ths)
356 {
357  if(ths->flags & PRECOMPUTE_WEIGHT)
358  nfft_free(ths->w);
359 
360  if(ths->flags & PRECOMPUTE_DAMP)
361  nfft_free(ths->w_hat);
362 
363  if(ths->flags & CGNR)
364  {
365  nfft_free(ths->v_iter);
366  nfft_free(ths->z_hat_iter);
367  }
368 
369  if(ths->flags & STEEPEST_DESCENT)
370  nfft_free(ths->v_iter);
371 
372  nfft_free(ths->p_hat_iter);
373  nfft_free(ths->f_hat_iter);
374 
375  nfft_free(ths->r_iter);
376  nfft_free(ths->y);
377 }
380 /****************************************************************************/
381 /****************************************************************************/
382 /****************************************************************************/
383 
384 void solver_init_advanced_double(solver_plan_double* ths, nfft_mv_plan_double *mv, unsigned flags)
385 {
386  ths->mv = mv;
387  ths->flags = flags;
388 
389  ths->y = (double*)nfft_malloc(ths->mv->M_total*sizeof(double));
390  ths->r_iter = (double*)nfft_malloc(ths->mv->M_total*sizeof(double));
391  ths->f_hat_iter = (double*)nfft_malloc(ths->mv->N_total*sizeof(double));
392  ths->p_hat_iter = (double*)nfft_malloc(ths->mv->N_total*sizeof(double));
393 
394  if(ths->flags & LANDWEBER)
395  ths->z_hat_iter = ths->p_hat_iter;
396 
397  if(ths->flags & STEEPEST_DESCENT)
398  {
399  ths->z_hat_iter = ths->p_hat_iter;
400  ths->v_iter = (double*)nfft_malloc(ths->mv->M_total*sizeof(double));
401  }
402 
403  if(ths->flags & CGNR)
404  {
405  ths->z_hat_iter = (double*)nfft_malloc(ths->mv->N_total*sizeof(double));
406  ths->v_iter = (double*)nfft_malloc(ths->mv->M_total*sizeof(double));
407  }
408 
409  if(ths->flags & CGNE)
410  ths->z_hat_iter = ths->p_hat_iter;
411 
412  if(ths->flags & PRECOMPUTE_WEIGHT)
413  ths->w = (double*) nfft_malloc(ths->mv->M_total*sizeof(double));
414 
415  if(ths->flags & PRECOMPUTE_DAMP)
416  ths->w_hat = (double*) nfft_malloc(ths->mv->N_total*sizeof(double));
417 }
418 
419 static void solver_init_double(solver_plan_double* ths, nfft_mv_plan_double *mv)
420 {
421  solver_init_advanced_double(ths, mv, CGNR);
422 }
423 
424 static void solver_before_loop_double(solver_plan_double* ths)
425 {
426  nfft_cp_double(ths->mv->f_hat, ths->f_hat_iter, ths->mv->N_total);
427 
428  NFFT_SWAP_double(ths->r_iter, ths->mv->f);
429  ths->mv->mv_trafo(ths->mv);
430  NFFT_SWAP_double(ths->r_iter, ths->mv->f);
431 
432  nfft_upd_axpy_double(ths->r_iter, -1.0, ths->y, ths->mv->M_total);
433 
434  if((!(ths->flags & LANDWEBER)) || (ths->flags & NORMS_FOR_LANDWEBER))
435  {
436  if(ths->flags & PRECOMPUTE_WEIGHT)
437  ths->dot_r_iter = nfft_dot_w_double(ths->r_iter, ths->w,
438  ths->mv->M_total);
439  else
440  ths->dot_r_iter = nfft_dot_double(ths->r_iter, ths->mv->M_total);
441  }
442 
443  /*-----------------*/
444  if(ths->flags & PRECOMPUTE_WEIGHT)
445  nfft_cp_w_double(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
446  else
447  nfft_cp_double(ths->mv->f, ths->r_iter, ths->mv->M_total);
448 
449  NFFT_SWAP_double(ths->z_hat_iter, ths->mv->f_hat);
450  ths->mv->mv_adjoint(ths->mv);
451  NFFT_SWAP_double(ths->z_hat_iter, ths->mv->f_hat);
452 
453  if((!(ths->flags & LANDWEBER)) || (ths->flags & NORMS_FOR_LANDWEBER))
454  {
455  if(ths->flags & PRECOMPUTE_DAMP)
456  ths->dot_z_hat_iter = nfft_dot_w_double(ths->z_hat_iter, ths->w_hat,
457  ths->mv->N_total);
458  else
459  ths->dot_z_hat_iter = nfft_dot_double(ths->z_hat_iter,
460  ths->mv->N_total);
461  }
462 
463  if(ths->flags & CGNE)
464  ths->dot_p_hat_iter = ths->dot_z_hat_iter;
465 
466  if(ths->flags & CGNR)
467  nfft_cp_double(ths->p_hat_iter, ths->z_hat_iter, ths->mv->N_total);
468 } /* void solver_before_loop */
469 
471 static void solver_loop_one_step_landweber_double(solver_plan_double* ths)
472 {
473  if(ths->flags & PRECOMPUTE_DAMP)
474  nfft_upd_xpawy_double(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
475  ths->z_hat_iter, ths->mv->N_total);
476  else
477  nfft_upd_xpay_double(ths->f_hat_iter, ths->alpha_iter, ths->z_hat_iter,
478  ths->mv->N_total);
479 
480  /*-----------------*/
481  nfft_cp_double(ths->mv->f_hat, ths->f_hat_iter, ths->mv->N_total);
482 
483  NFFT_SWAP_double(ths->r_iter,ths->mv->f);
484  ths->mv->mv_trafo(ths->mv);
485  NFFT_SWAP_double(ths->r_iter,ths->mv->f);
486 
487  nfft_upd_axpy_double(ths->r_iter, -1.0, ths->y, ths->mv->M_total);
488 
489  if(ths->flags & NORMS_FOR_LANDWEBER)
490  {
491  if(ths->flags & PRECOMPUTE_WEIGHT)
492  ths->dot_r_iter = nfft_dot_w_double(ths->r_iter,ths->w,
493  ths->mv->M_total);
494  else
495  ths->dot_r_iter = nfft_dot_double(ths->r_iter, ths->mv->M_total);
496  }
497 
498  /*-----------------*/
499  if(ths->flags & PRECOMPUTE_WEIGHT)
500  nfft_cp_w_double(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
501  else
502  nfft_cp_double(ths->mv->f, ths->r_iter, ths->mv->M_total);
503 
504  NFFT_SWAP_double(ths->z_hat_iter,ths->mv->f_hat);
505  ths->mv->mv_adjoint(ths->mv);
506  NFFT_SWAP_double(ths->z_hat_iter,ths->mv->f_hat);
507 
508  if(ths->flags & NORMS_FOR_LANDWEBER)
509  {
510  if(ths->flags & PRECOMPUTE_DAMP)
511  ths->dot_z_hat_iter = nfft_dot_w_double(ths->z_hat_iter, ths->w_hat,
512  ths->mv->N_total);
513  else
514  ths->dot_z_hat_iter = nfft_dot_double(ths->z_hat_iter,
515  ths->mv->N_total);
516  }
517 } /* void solver_loop_one_step_landweber */
518 
520 static void solver_loop_one_step_steepest_descent_double(solver_plan_double *ths)
521 {
522  if(ths->flags & PRECOMPUTE_DAMP)
523  nfft_cp_w_double(ths->mv->f_hat, ths->w_hat, ths->z_hat_iter,
524  ths->mv->N_total);
525  else
526  nfft_cp_double(ths->mv->f_hat, ths->z_hat_iter, ths->mv->N_total);
527 
528  NFFT_SWAP_double(ths->v_iter,ths->mv->f);
529  ths->mv->mv_trafo(ths->mv);
530  NFFT_SWAP_double(ths->v_iter,ths->mv->f);
531 
532  if(ths->flags & PRECOMPUTE_WEIGHT)
533  ths->dot_v_iter = nfft_dot_w_double(ths->v_iter,ths->w,ths->mv->M_total);
534  else
535  ths->dot_v_iter = nfft_dot_double(ths->v_iter, ths->mv->M_total);
536 
537  /*-----------------*/
538  ths->alpha_iter = ths->dot_z_hat_iter / ths->dot_v_iter;
539 
540  /*-----------------*/
541  if(ths->flags & PRECOMPUTE_DAMP)
542  nfft_upd_xpawy_double(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
543  ths->z_hat_iter, ths->mv->N_total);
544  else
545  nfft_upd_xpay_double(ths->f_hat_iter, ths->alpha_iter, ths->z_hat_iter,
546  ths->mv->N_total);
547 
548  /*-----------------*/
549  nfft_upd_xpay_double(ths->r_iter, -ths->alpha_iter, ths->v_iter,
550  ths->mv->M_total);
551 
552  if(ths->flags & PRECOMPUTE_WEIGHT)
553  ths->dot_r_iter = nfft_dot_w_double(ths->r_iter,ths->w,ths->mv->M_total);
554  else
555  ths->dot_r_iter = nfft_dot_double(ths->r_iter, ths->mv->M_total);
556 
557  /*-----------------*/
558  if(ths->flags & PRECOMPUTE_WEIGHT)
559  nfft_cp_w_double(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
560  else
561  nfft_cp_double(ths->mv->f, ths->r_iter, ths->mv->M_total);
562 
563  NFFT_SWAP_double(ths->z_hat_iter,ths->mv->f_hat);
564  ths->mv->mv_adjoint(ths->mv);
565  NFFT_SWAP_double(ths->z_hat_iter,ths->mv->f_hat);
566 
567  if(ths->flags & PRECOMPUTE_DAMP)
568  ths->dot_z_hat_iter = nfft_dot_w_double(ths->z_hat_iter, ths->w_hat,
569  ths->mv->N_total);
570  else
571  ths->dot_z_hat_iter = nfft_dot_double(ths->z_hat_iter, ths->mv->N_total);
572 } /* void solver_loop_one_step_steepest_descent */
573 
575 static void solver_loop_one_step_cgnr_double(solver_plan_double *ths)
576 {
577  if(ths->flags & PRECOMPUTE_DAMP)
578  nfft_cp_w_double(ths->mv->f_hat, ths->w_hat, ths->p_hat_iter,
579  ths->mv->N_total);
580  else
581  nfft_cp_double(ths->mv->f_hat, ths->p_hat_iter, ths->mv->N_total);
582 
583  NFFT_SWAP_double(ths->v_iter,ths->mv->f);
584  ths->mv->mv_trafo(ths->mv);
585  NFFT_SWAP_double(ths->v_iter,ths->mv->f);
586 
587  if(ths->flags & PRECOMPUTE_WEIGHT)
588  ths->dot_v_iter = nfft_dot_w_double(ths->v_iter,ths->w,ths->mv->M_total);
589  else
590  ths->dot_v_iter = nfft_dot_double(ths->v_iter, ths->mv->M_total);
591 
592  /*-----------------*/
593  ths->alpha_iter = ths->dot_z_hat_iter / ths->dot_v_iter;
594 
595  /*-----------------*/
596  if(ths->flags & PRECOMPUTE_DAMP)
597  nfft_upd_xpawy_double(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
598  ths->p_hat_iter, ths->mv->N_total);
599  else
600  nfft_upd_xpay_double(ths->f_hat_iter, ths->alpha_iter, ths->p_hat_iter,
601  ths->mv->N_total);
602 
603  /*-----------------*/
604  nfft_upd_xpay_double(ths->r_iter, -ths->alpha_iter, ths->v_iter,
605  ths->mv->M_total);
606 
607  if(ths->flags & PRECOMPUTE_WEIGHT)
608  ths->dot_r_iter = nfft_dot_w_double(ths->r_iter,ths->w,ths->mv->M_total);
609  else
610  ths->dot_r_iter = nfft_dot_double(ths->r_iter, ths->mv->M_total);
611 
612  /*-----------------*/
613  if(ths->flags & PRECOMPUTE_WEIGHT)
614  nfft_cp_w_double(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
615  else
616  nfft_cp_double(ths->mv->f, ths->r_iter, ths->mv->M_total);
617 
618  NFFT_SWAP_double(ths->z_hat_iter,ths->mv->f_hat);
619  ths->mv->mv_adjoint(ths->mv);
620  NFFT_SWAP_double(ths->z_hat_iter,ths->mv->f_hat);
621 
622  ths->dot_z_hat_iter_old = ths->dot_z_hat_iter;
623  if(ths->flags & PRECOMPUTE_DAMP)
624  ths->dot_z_hat_iter = nfft_dot_w_double(ths->z_hat_iter, ths->w_hat,
625  ths->mv->N_total);
626  else
627  ths->dot_z_hat_iter = nfft_dot_double(ths->z_hat_iter, ths->mv->N_total);
628 
629  /*-----------------*/
630  ths->beta_iter = ths->dot_z_hat_iter / ths->dot_z_hat_iter_old;
631 
632  /*-----------------*/
633  nfft_upd_axpy_double(ths->p_hat_iter, ths->beta_iter, ths->z_hat_iter,
634  ths->mv->N_total);
635 } /* void solver_loop_one_step_cgnr */
636 
638 static void solver_loop_one_step_cgne_double(solver_plan_double *ths)
639 {
640  ths->alpha_iter = ths->dot_r_iter / ths->dot_p_hat_iter;
641 
642  /*-----------------*/
643  if(ths->flags & PRECOMPUTE_DAMP)
644  nfft_upd_xpawy_double(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
645  ths->p_hat_iter, ths->mv->N_total);
646  else
647  nfft_upd_xpay_double(ths->f_hat_iter, ths->alpha_iter, ths->p_hat_iter,
648  ths->mv->N_total);
649 
650  /*-----------------*/
651  if(ths->flags & PRECOMPUTE_DAMP)
652  nfft_cp_w_double(ths->mv->f_hat, ths->w_hat, ths->p_hat_iter,
653  ths->mv->N_total);
654  else
655  nfft_cp_double(ths->mv->f_hat, ths->p_hat_iter, ths->mv->N_total);
656 
657  ths->mv->mv_trafo(ths->mv);
658 
659  nfft_upd_xpay_double(ths->r_iter, -ths->alpha_iter, ths->mv->f,
660  ths->mv->M_total);
661 
662  ths->dot_r_iter_old = ths->dot_r_iter;
663  if(ths->flags & PRECOMPUTE_WEIGHT)
664  ths->dot_r_iter = nfft_dot_w_double(ths->r_iter,ths->w,ths->mv->M_total);
665  else
666  ths->dot_r_iter = nfft_dot_double(ths->r_iter, ths->mv->M_total);
667 
668  /*-----------------*/
669  ths->beta_iter = ths->dot_r_iter / ths->dot_r_iter_old;
670 
671  /*-----------------*/
672  if(ths->flags & PRECOMPUTE_WEIGHT)
673  nfft_cp_w_double(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
674  else
675  nfft_cp_double(ths->mv->f, ths->r_iter, ths->mv->M_total);
676 
677  ths->mv->mv_adjoint(ths->mv);
678 
679  nfft_upd_axpy_double(ths->p_hat_iter, ths->beta_iter, ths->mv->f_hat,
680  ths->mv->N_total);
681 
682  if(ths->flags & PRECOMPUTE_DAMP)
683  ths->dot_p_hat_iter = nfft_dot_w_double(ths->p_hat_iter, ths->w_hat,
684  ths->mv->N_total);
685  else
686  ths->dot_p_hat_iter = nfft_dot_double(ths->p_hat_iter, ths->mv->N_total);
687 } /* void solver_loop_one_step_cgne */
688 
690 static void solver_loop_one_step_double(solver_plan_double *ths)
691 {
692  if(ths->flags & LANDWEBER)
693  solver_loop_one_step_landweber_double(ths);
694 
695  if(ths->flags & STEEPEST_DESCENT)
696  solver_loop_one_step_steepest_descent_double(ths);
697 
698  if(ths->flags & CGNR)
699  solver_loop_one_step_cgnr_double(ths);
700 
701  if(ths->flags & CGNE)
702  solver_loop_one_step_cgne_double(ths);
703 } /* void solver_loop_one_step */
704 
706 static void solver_finalize_double(solver_plan_double *ths)
707 {
708  if(ths->flags & PRECOMPUTE_WEIGHT)
709  nfft_free(ths->w);
710 
711  if(ths->flags & PRECOMPUTE_DAMP)
712  nfft_free(ths->w_hat);
713 
714  if(ths->flags & CGNR)
715  {
716  nfft_free(ths->v_iter);
717  nfft_free(ths->z_hat_iter);
718  }
719 
720  if(ths->flags & STEEPEST_DESCENT)
721  nfft_free(ths->v_iter);
722 
723  nfft_free(ths->p_hat_iter);
724  nfft_free(ths->f_hat_iter);
725 
726  nfft_free(ths->r_iter);
727  nfft_free(ths->y);
728 }
Generated on Tue Apr 30 2013 by Doxygen 1.8.1