NFFT  3.4.1
solver.c
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2002, 2017 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 
23 #include "config.h"
24 
25 #ifdef HAVE_COMPLEX_H
26 #include <complex.h>
27 #endif
28 #include "nfft3.h"
29 #include "infft.h"
30 
31 #undef X
32 #if defined(NFFT_SINGLE)
33 #define X(name) CONCAT(solverf_,name)
34 #elif defined(NFFT_LDOUBLE)
35 #define X(name) CONCAT(solverl_,name)
36 #else
37 #define X(name) CONCAT(solver_,name)
38 #endif
39 
40 void X(init_advanced_complex)(X(plan_complex)* ths, Y(mv_plan_complex) *mv,
41  unsigned flags)
42 {
43  ths->mv = mv;
44  ths->flags = flags;
45 
46  ths->y = (C*)Y(malloc)((size_t)(ths->mv->M_total) * sizeof(C));
47  ths->r_iter = (C*)Y(malloc)((size_t)(ths->mv->M_total) * sizeof(C));
48  ths->f_hat_iter = (C*)Y(malloc)((size_t)(ths->mv->N_total) * sizeof(C));
49  ths->p_hat_iter = (C*)Y(malloc)((size_t)(ths->mv->N_total) * sizeof(C));
50 
51  if(ths->flags & LANDWEBER)
52  ths->z_hat_iter = ths->p_hat_iter;
53 
54  if(ths->flags & STEEPEST_DESCENT)
55  {
56  ths->z_hat_iter = ths->p_hat_iter;
57  ths->v_iter = (C*)Y(malloc)((size_t)(ths->mv->M_total) * sizeof(C));
58  }
59 
60  if(ths->flags & CGNR)
61  {
62  ths->z_hat_iter = (C*)Y(malloc)((size_t)(ths->mv->N_total) * sizeof(C));
63  ths->v_iter = (C*)Y(malloc)((size_t)(ths->mv->M_total) * sizeof(C));
64  }
65 
66  if(ths->flags & CGNE)
67  ths->z_hat_iter = ths->p_hat_iter;
68 
69  if(ths->flags & PRECOMPUTE_WEIGHT)
70  ths->w = (R*) Y(malloc)((size_t)(ths->mv->M_total) * sizeof(R));
71 
72  if(ths->flags & PRECOMPUTE_DAMP)
73  ths->w_hat = (R*) Y(malloc)((size_t)(ths->mv->N_total) * sizeof(R));
74 }
75 
76 void X(init_complex)(X(plan_complex)* ths, Y(mv_plan_complex) *mv)
77 {
78  X(init_advanced_complex)(ths, mv, CGNR);
79 }
80 
81 void X(before_loop_complex)(X(plan_complex)* ths)
82 {
83  Y(cp_complex)(ths->mv->f_hat, ths->f_hat_iter, ths->mv->N_total);
84 
85  CSWAP(ths->r_iter, ths->mv->f);
86  ths->mv->mv_trafo(ths->mv);
87  CSWAP(ths->r_iter, ths->mv->f);
88 
89  Y(upd_axpy_complex)(ths->r_iter, K(-1.0), ths->y, ths->mv->M_total);
90 
91  if((!(ths->flags & LANDWEBER)) || (ths->flags & NORMS_FOR_LANDWEBER))
92  {
93  if(ths->flags & PRECOMPUTE_WEIGHT)
94  ths->dot_r_iter = Y(dot_w_complex)(ths->r_iter, ths->w,
95  ths->mv->M_total);
96  else
97  ths->dot_r_iter = Y(dot_complex)(ths->r_iter, ths->mv->M_total);
98  }
99 
100  /*-----------------*/
101  if(ths->flags & PRECOMPUTE_WEIGHT)
102  Y(cp_w_complex)(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
103  else
104  Y(cp_complex)(ths->mv->f, ths->r_iter, ths->mv->M_total);
105 
106  CSWAP(ths->z_hat_iter, ths->mv->f_hat);
107  ths->mv->mv_adjoint(ths->mv);
108  CSWAP(ths->z_hat_iter, ths->mv->f_hat);
109 
110  if((!(ths->flags & LANDWEBER)) || (ths->flags & NORMS_FOR_LANDWEBER))
111  {
112  if(ths->flags & PRECOMPUTE_DAMP)
113  ths->dot_z_hat_iter = Y(dot_w_complex)(ths->z_hat_iter, ths->w_hat,
114  ths->mv->N_total);
115  else
116  ths->dot_z_hat_iter = Y(dot_complex)(ths->z_hat_iter,
117  ths->mv->N_total);
118  }
119 
120  if(ths->flags & CGNE)
121  ths->dot_p_hat_iter = ths->dot_z_hat_iter;
122 
123  if(ths->flags & CGNR)
124  Y(cp_complex)(ths->p_hat_iter, ths->z_hat_iter, ths->mv->N_total);
125 } /* void solver_before_loop */
126 
128 static void solver_loop_one_step_landweber_complex(X(plan_complex)* ths)
129 {
130  if(ths->flags & PRECOMPUTE_DAMP)
131  Y(upd_xpawy_complex)(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
132  ths->z_hat_iter, ths->mv->N_total);
133  else
134  Y(upd_xpay_complex)(ths->f_hat_iter, ths->alpha_iter, ths->z_hat_iter,
135  ths->mv->N_total);
136 
137  /*-----------------*/
138  Y(cp_complex)(ths->mv->f_hat, ths->f_hat_iter, ths->mv->N_total);
139 
140  CSWAP(ths->r_iter,ths->mv->f);
141  ths->mv->mv_trafo(ths->mv);
142  CSWAP(ths->r_iter,ths->mv->f);
143 
144  Y(upd_axpy_complex)(ths->r_iter, K(-1.0), ths->y, ths->mv->M_total);
145 
146  if(ths->flags & NORMS_FOR_LANDWEBER)
147  {
148  if(ths->flags & PRECOMPUTE_WEIGHT)
149  ths->dot_r_iter = Y(dot_w_complex)(ths->r_iter,ths->w,
150  ths->mv->M_total);
151  else
152  ths->dot_r_iter = Y(dot_complex)(ths->r_iter, ths->mv->M_total);
153  }
154 
155  /*-----------------*/
156  if(ths->flags & PRECOMPUTE_WEIGHT)
157  Y(cp_w_complex)(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
158  else
159  Y(cp_complex)(ths->mv->f, ths->r_iter, ths->mv->M_total);
160 
161  CSWAP(ths->z_hat_iter,ths->mv->f_hat);
162  ths->mv->mv_adjoint(ths->mv);
163  CSWAP(ths->z_hat_iter,ths->mv->f_hat);
164 
165  if(ths->flags & NORMS_FOR_LANDWEBER)
166  {
167  if(ths->flags & PRECOMPUTE_DAMP)
168  ths->dot_z_hat_iter = Y(dot_w_complex)(ths->z_hat_iter, ths->w_hat,
169  ths->mv->N_total);
170  else
171  ths->dot_z_hat_iter = Y(dot_complex)(ths->z_hat_iter,
172  ths->mv->N_total);
173  }
174 } /* void solver_loop_one_step_landweber */
175 
177 static void solver_loop_one_step_steepest_descent_complex(X(plan_complex) *ths)
178 {
179  if(ths->flags & PRECOMPUTE_DAMP)
180  Y(cp_w_complex)(ths->mv->f_hat, ths->w_hat, ths->z_hat_iter,
181  ths->mv->N_total);
182  else
183  Y(cp_complex)(ths->mv->f_hat, ths->z_hat_iter, ths->mv->N_total);
184 
185  CSWAP(ths->v_iter,ths->mv->f);
186  ths->mv->mv_trafo(ths->mv);
187  CSWAP(ths->v_iter,ths->mv->f);
188 
189  if(ths->flags & PRECOMPUTE_WEIGHT)
190  ths->dot_v_iter = Y(dot_w_complex)(ths->v_iter,ths->w,ths->mv->M_total);
191  else
192  ths->dot_v_iter = Y(dot_complex)(ths->v_iter, ths->mv->M_total);
193 
194  /*-----------------*/
195  ths->alpha_iter = ths->dot_z_hat_iter / ths->dot_v_iter;
196 
197  /*-----------------*/
198  if(ths->flags & PRECOMPUTE_DAMP)
199  Y(upd_xpawy_complex)(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
200  ths->z_hat_iter, ths->mv->N_total);
201  else
202  Y(upd_xpay_complex)(ths->f_hat_iter, ths->alpha_iter, ths->z_hat_iter,
203  ths->mv->N_total);
204 
205  /*-----------------*/
206  Y(upd_xpay_complex)(ths->r_iter, -ths->alpha_iter, ths->v_iter,
207  ths->mv->M_total);
208 
209  if(ths->flags & PRECOMPUTE_WEIGHT)
210  ths->dot_r_iter = Y(dot_w_complex)(ths->r_iter,ths->w,ths->mv->M_total);
211  else
212  ths->dot_r_iter = Y(dot_complex)(ths->r_iter, ths->mv->M_total);
213 
214  /*-----------------*/
215  if(ths->flags & PRECOMPUTE_WEIGHT)
216  Y(cp_w_complex)(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
217  else
218  Y(cp_complex)(ths->mv->f, ths->r_iter, ths->mv->M_total);
219 
220  CSWAP(ths->z_hat_iter,ths->mv->f_hat);
221  ths->mv->mv_adjoint(ths->mv);
222  CSWAP(ths->z_hat_iter,ths->mv->f_hat);
223 
224  if(ths->flags & PRECOMPUTE_DAMP)
225  ths->dot_z_hat_iter = Y(dot_w_complex)(ths->z_hat_iter, ths->w_hat,
226  ths->mv->N_total);
227  else
228  ths->dot_z_hat_iter = Y(dot_complex)(ths->z_hat_iter, ths->mv->N_total);
229 } /* void solver_loop_one_step_steepest_descent */
230 
232 static void solver_loop_one_step_cgnr_complex(X(plan_complex) *ths)
233 {
234  if(ths->flags & PRECOMPUTE_DAMP)
235  Y(cp_w_complex)(ths->mv->f_hat, ths->w_hat, ths->p_hat_iter,
236  ths->mv->N_total);
237  else
238  Y(cp_complex)(ths->mv->f_hat, ths->p_hat_iter, ths->mv->N_total);
239 
240  CSWAP(ths->v_iter,ths->mv->f);
241  ths->mv->mv_trafo(ths->mv);
242  CSWAP(ths->v_iter,ths->mv->f);
243 
244  if(ths->flags & PRECOMPUTE_WEIGHT)
245  ths->dot_v_iter = Y(dot_w_complex)(ths->v_iter,ths->w,ths->mv->M_total);
246  else
247  ths->dot_v_iter = Y(dot_complex)(ths->v_iter, ths->mv->M_total);
248 
249  /*-----------------*/
250  ths->alpha_iter = ths->dot_z_hat_iter / ths->dot_v_iter;
251 
252  /*-----------------*/
253  if(ths->flags & PRECOMPUTE_DAMP)
254  Y(upd_xpawy_complex)(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
255  ths->p_hat_iter, ths->mv->N_total);
256  else
257  Y(upd_xpay_complex)(ths->f_hat_iter, ths->alpha_iter, ths->p_hat_iter,
258  ths->mv->N_total);
259 
260  /*-----------------*/
261  Y(upd_xpay_complex)(ths->r_iter, -ths->alpha_iter, ths->v_iter,
262  ths->mv->M_total);
263 
264  if(ths->flags & PRECOMPUTE_WEIGHT)
265  ths->dot_r_iter = Y(dot_w_complex)(ths->r_iter,ths->w,ths->mv->M_total);
266  else
267  ths->dot_r_iter = Y(dot_complex)(ths->r_iter, ths->mv->M_total);
268 
269  /*-----------------*/
270  if(ths->flags & PRECOMPUTE_WEIGHT)
271  Y(cp_w_complex)(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
272  else
273  Y(cp_complex)(ths->mv->f, ths->r_iter, ths->mv->M_total);
274 
275  CSWAP(ths->z_hat_iter,ths->mv->f_hat);
276  ths->mv->mv_adjoint(ths->mv);
277  CSWAP(ths->z_hat_iter,ths->mv->f_hat);
278 
279  ths->dot_z_hat_iter_old = ths->dot_z_hat_iter;
280  if(ths->flags & PRECOMPUTE_DAMP)
281  ths->dot_z_hat_iter = Y(dot_w_complex)(ths->z_hat_iter, ths->w_hat,
282  ths->mv->N_total);
283  else
284  ths->dot_z_hat_iter = Y(dot_complex)(ths->z_hat_iter, ths->mv->N_total);
285 
286  /*-----------------*/
287  ths->beta_iter = ths->dot_z_hat_iter / ths->dot_z_hat_iter_old;
288 
289  /*-----------------*/
290  Y(upd_axpy_complex)(ths->p_hat_iter, ths->beta_iter, ths->z_hat_iter,
291  ths->mv->N_total);
292 } /* void solver_loop_one_step_cgnr */
293 
295 static void solver_loop_one_step_cgne_complex(X(plan_complex) *ths)
296 {
297  ths->alpha_iter = ths->dot_r_iter / ths->dot_p_hat_iter;
298 
299  /*-----------------*/
300  if(ths->flags & PRECOMPUTE_DAMP)
301  Y(upd_xpawy_complex)(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
302  ths->p_hat_iter, ths->mv->N_total);
303  else
304  Y(upd_xpay_complex)(ths->f_hat_iter, ths->alpha_iter, ths->p_hat_iter,
305  ths->mv->N_total);
306 
307  /*-----------------*/
308  if(ths->flags & PRECOMPUTE_DAMP)
309  Y(cp_w_complex)(ths->mv->f_hat, ths->w_hat, ths->p_hat_iter,
310  ths->mv->N_total);
311  else
312  Y(cp_complex)(ths->mv->f_hat, ths->p_hat_iter, ths->mv->N_total);
313 
314  ths->mv->mv_trafo(ths->mv);
315 
316  Y(upd_xpay_complex)(ths->r_iter, -ths->alpha_iter, ths->mv->f,
317  ths->mv->M_total);
318 
319  ths->dot_r_iter_old = ths->dot_r_iter;
320  if(ths->flags & PRECOMPUTE_WEIGHT)
321  ths->dot_r_iter = Y(dot_w_complex)(ths->r_iter,ths->w,ths->mv->M_total);
322  else
323  ths->dot_r_iter = Y(dot_complex)(ths->r_iter, ths->mv->M_total);
324 
325  /*-----------------*/
326  ths->beta_iter = ths->dot_r_iter / ths->dot_r_iter_old;
327 
328  /*-----------------*/
329  if(ths->flags & PRECOMPUTE_WEIGHT)
330  Y(cp_w_complex)(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
331  else
332  Y(cp_complex)(ths->mv->f, ths->r_iter, ths->mv->M_total);
333 
334  ths->mv->mv_adjoint(ths->mv);
335 
336  Y(upd_axpy_complex)(ths->p_hat_iter, ths->beta_iter, ths->mv->f_hat,
337  ths->mv->N_total);
338 
339  if(ths->flags & PRECOMPUTE_DAMP)
340  ths->dot_p_hat_iter = Y(dot_w_complex)(ths->p_hat_iter, ths->w_hat,
341  ths->mv->N_total);
342  else
343  ths->dot_p_hat_iter = Y(dot_complex)(ths->p_hat_iter, ths->mv->N_total);
344 } /* void solver_loop_one_step_cgne */
345 
347 void X(loop_one_step_complex)(X(plan_complex) *ths)
348 {
349  if(ths->flags & LANDWEBER)
350  solver_loop_one_step_landweber_complex(ths);
351 
352  if(ths->flags & STEEPEST_DESCENT)
353  solver_loop_one_step_steepest_descent_complex(ths);
354 
355  if(ths->flags & CGNR)
356  solver_loop_one_step_cgnr_complex(ths);
357 
358  if(ths->flags & CGNE)
359  solver_loop_one_step_cgne_complex(ths);
360 } /* void solver_loop_one_step */
361 
363 void X(finalize_complex)(X(plan_complex) *ths)
364 {
365  if(ths->flags & PRECOMPUTE_WEIGHT)
366  Y(free)(ths->w);
367 
368  if(ths->flags & PRECOMPUTE_DAMP)
369  Y(free)(ths->w_hat);
370 
371  if(ths->flags & CGNR)
372  {
373  Y(free)(ths->v_iter);
374  Y(free)(ths->z_hat_iter);
375  }
376 
377  if(ths->flags & STEEPEST_DESCENT)
378  Y(free)(ths->v_iter);
379 
380  Y(free)(ths->p_hat_iter);
381  Y(free)(ths->f_hat_iter);
382 
383  Y(free)(ths->r_iter);
384  Y(free)(ths->y);
385 }
388 /****************************************************************************/
389 /****************************************************************************/
390 /****************************************************************************/
391 
392 void X(init_advanced_double)(X(plan_double)* ths, Y(mv_plan_double) *mv, unsigned flags)
393 {
394  ths->mv = mv;
395  ths->flags = flags;
396 
397  ths->y = (R*)Y(malloc)((size_t)(ths->mv->M_total) * sizeof(R));
398  ths->r_iter = (R*)Y(malloc)((size_t)(ths->mv->M_total) * sizeof(R));
399  ths->f_hat_iter = (R*)Y(malloc)((size_t)(ths->mv->N_total) * sizeof(R));
400  ths->p_hat_iter = (R*)Y(malloc)((size_t)(ths->mv->N_total) * sizeof(R));
401 
402  if(ths->flags & LANDWEBER)
403  ths->z_hat_iter = ths->p_hat_iter;
404 
405  if(ths->flags & STEEPEST_DESCENT)
406  {
407  ths->z_hat_iter = ths->p_hat_iter;
408  ths->v_iter = (R*)Y(malloc)((size_t)(ths->mv->M_total) * sizeof(R));
409  }
410 
411  if(ths->flags & CGNR)
412  {
413  ths->z_hat_iter = (R*)Y(malloc)((size_t)(ths->mv->N_total) * sizeof(R));
414  ths->v_iter = (R*)Y(malloc)((size_t)(ths->mv->M_total) * sizeof(R));
415  }
416 
417  if(ths->flags & CGNE)
418  ths->z_hat_iter = ths->p_hat_iter;
419 
420  if(ths->flags & PRECOMPUTE_WEIGHT)
421  ths->w = (R*) Y(malloc)((size_t)(ths->mv->M_total) * sizeof(R));
422 
423  if(ths->flags & PRECOMPUTE_DAMP)
424  ths->w_hat = (R*) Y(malloc)((size_t)(ths->mv->N_total) * sizeof(R));
425 }
426 
427 void X(init_double)(X(plan_double)* ths, Y(mv_plan_double) *mv)
428 {
429  X(init_advanced_double)(ths, mv, CGNR);
430 }
431 
432 void X(before_loop_double)(X(plan_double)* ths)
433 {
434  Y(cp_double)(ths->mv->f_hat, ths->f_hat_iter, ths->mv->N_total);
435 
436  RSWAP(ths->r_iter, ths->mv->f);
437  ths->mv->mv_trafo(ths->mv);
438  RSWAP(ths->r_iter, ths->mv->f);
439 
440  Y(upd_axpy_double)(ths->r_iter, K(-1.0), ths->y, ths->mv->M_total);
441 
442  if((!(ths->flags & LANDWEBER)) || (ths->flags & NORMS_FOR_LANDWEBER))
443  {
444  if(ths->flags & PRECOMPUTE_WEIGHT)
445  ths->dot_r_iter = Y(dot_w_double)(ths->r_iter, ths->w,
446  ths->mv->M_total);
447  else
448  ths->dot_r_iter = Y(dot_double)(ths->r_iter, ths->mv->M_total);
449  }
450 
451  /*-----------------*/
452  if(ths->flags & PRECOMPUTE_WEIGHT)
453  Y(cp_w_double)(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
454  else
455  Y(cp_double)(ths->mv->f, ths->r_iter, ths->mv->M_total);
456 
457  RSWAP(ths->z_hat_iter, ths->mv->f_hat);
458  ths->mv->mv_adjoint(ths->mv);
459  RSWAP(ths->z_hat_iter, ths->mv->f_hat);
460 
461  if((!(ths->flags & LANDWEBER)) || (ths->flags & NORMS_FOR_LANDWEBER))
462  {
463  if(ths->flags & PRECOMPUTE_DAMP)
464  ths->dot_z_hat_iter = Y(dot_w_double)(ths->z_hat_iter, ths->w_hat,
465  ths->mv->N_total);
466  else
467  ths->dot_z_hat_iter = Y(dot_double)(ths->z_hat_iter,
468  ths->mv->N_total);
469  }
470 
471  if(ths->flags & CGNE)
472  ths->dot_p_hat_iter = ths->dot_z_hat_iter;
473 
474  if(ths->flags & CGNR)
475  Y(cp_double)(ths->p_hat_iter, ths->z_hat_iter, ths->mv->N_total);
476 } /* void solver_before_loop */
477 
479 static void solver_loop_one_step_landweber_double(X(plan_double)* ths)
480 {
481  if(ths->flags & PRECOMPUTE_DAMP)
482  Y(upd_xpawy_double)(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
483  ths->z_hat_iter, ths->mv->N_total);
484  else
485  Y(upd_xpay_double)(ths->f_hat_iter, ths->alpha_iter, ths->z_hat_iter,
486  ths->mv->N_total);
487 
488  /*-----------------*/
489  Y(cp_double)(ths->mv->f_hat, ths->f_hat_iter, ths->mv->N_total);
490 
491  RSWAP(ths->r_iter,ths->mv->f);
492  ths->mv->mv_trafo(ths->mv);
493  RSWAP(ths->r_iter,ths->mv->f);
494 
495  Y(upd_axpy_double)(ths->r_iter, K(-1.0), ths->y, ths->mv->M_total);
496 
497  if(ths->flags & NORMS_FOR_LANDWEBER)
498  {
499  if(ths->flags & PRECOMPUTE_WEIGHT)
500  ths->dot_r_iter = Y(dot_w_double)(ths->r_iter,ths->w,
501  ths->mv->M_total);
502  else
503  ths->dot_r_iter = Y(dot_double)(ths->r_iter, ths->mv->M_total);
504  }
505 
506  /*-----------------*/
507  if(ths->flags & PRECOMPUTE_WEIGHT)
508  Y(cp_w_double)(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
509  else
510  Y(cp_double)(ths->mv->f, ths->r_iter, ths->mv->M_total);
511 
512  RSWAP(ths->z_hat_iter,ths->mv->f_hat);
513  ths->mv->mv_adjoint(ths->mv);
514  RSWAP(ths->z_hat_iter,ths->mv->f_hat);
515 
516  if(ths->flags & NORMS_FOR_LANDWEBER)
517  {
518  if(ths->flags & PRECOMPUTE_DAMP)
519  ths->dot_z_hat_iter = Y(dot_w_double)(ths->z_hat_iter, ths->w_hat,
520  ths->mv->N_total);
521  else
522  ths->dot_z_hat_iter = Y(dot_double)(ths->z_hat_iter,
523  ths->mv->N_total);
524  }
525 } /* void solver_loop_one_step_landweber */
526 
528 static void solver_loop_one_step_steepest_descent_double(X(plan_double) *ths)
529 {
530  if(ths->flags & PRECOMPUTE_DAMP)
531  Y(cp_w_double)(ths->mv->f_hat, ths->w_hat, ths->z_hat_iter,
532  ths->mv->N_total);
533  else
534  Y(cp_double)(ths->mv->f_hat, ths->z_hat_iter, ths->mv->N_total);
535 
536  RSWAP(ths->v_iter,ths->mv->f);
537  ths->mv->mv_trafo(ths->mv);
538  RSWAP(ths->v_iter,ths->mv->f);
539 
540  if(ths->flags & PRECOMPUTE_WEIGHT)
541  ths->dot_v_iter = Y(dot_w_double)(ths->v_iter,ths->w,ths->mv->M_total);
542  else
543  ths->dot_v_iter = Y(dot_double)(ths->v_iter, ths->mv->M_total);
544 
545  /*-----------------*/
546  ths->alpha_iter = ths->dot_z_hat_iter / ths->dot_v_iter;
547 
548  /*-----------------*/
549  if(ths->flags & PRECOMPUTE_DAMP)
550  Y(upd_xpawy_double)(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
551  ths->z_hat_iter, ths->mv->N_total);
552  else
553  Y(upd_xpay_double)(ths->f_hat_iter, ths->alpha_iter, ths->z_hat_iter,
554  ths->mv->N_total);
555 
556  /*-----------------*/
557  Y(upd_xpay_double)(ths->r_iter, -ths->alpha_iter, ths->v_iter,
558  ths->mv->M_total);
559 
560  if(ths->flags & PRECOMPUTE_WEIGHT)
561  ths->dot_r_iter = Y(dot_w_double)(ths->r_iter,ths->w,ths->mv->M_total);
562  else
563  ths->dot_r_iter = Y(dot_double)(ths->r_iter, ths->mv->M_total);
564 
565  /*-----------------*/
566  if(ths->flags & PRECOMPUTE_WEIGHT)
567  Y(cp_w_double)(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
568  else
569  Y(cp_double)(ths->mv->f, ths->r_iter, ths->mv->M_total);
570 
571  RSWAP(ths->z_hat_iter,ths->mv->f_hat);
572  ths->mv->mv_adjoint(ths->mv);
573  RSWAP(ths->z_hat_iter,ths->mv->f_hat);
574 
575  if(ths->flags & PRECOMPUTE_DAMP)
576  ths->dot_z_hat_iter = Y(dot_w_double)(ths->z_hat_iter, ths->w_hat,
577  ths->mv->N_total);
578  else
579  ths->dot_z_hat_iter = Y(dot_double)(ths->z_hat_iter, ths->mv->N_total);
580 } /* void solver_loop_one_step_steepest_descent */
581 
583 static void solver_loop_one_step_cgnr_double(X(plan_double) *ths)
584 {
585  if(ths->flags & PRECOMPUTE_DAMP)
586  Y(cp_w_double)(ths->mv->f_hat, ths->w_hat, ths->p_hat_iter,
587  ths->mv->N_total);
588  else
589  Y(cp_double)(ths->mv->f_hat, ths->p_hat_iter, ths->mv->N_total);
590 
591  RSWAP(ths->v_iter,ths->mv->f);
592  ths->mv->mv_trafo(ths->mv);
593  RSWAP(ths->v_iter,ths->mv->f);
594 
595  if(ths->flags & PRECOMPUTE_WEIGHT)
596  ths->dot_v_iter = Y(dot_w_double)(ths->v_iter,ths->w,ths->mv->M_total);
597  else
598  ths->dot_v_iter = Y(dot_double)(ths->v_iter, ths->mv->M_total);
599 
600  /*-----------------*/
601  ths->alpha_iter = ths->dot_z_hat_iter / ths->dot_v_iter;
602 
603  /*-----------------*/
604  if(ths->flags & PRECOMPUTE_DAMP)
605  Y(upd_xpawy_double)(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
606  ths->p_hat_iter, ths->mv->N_total);
607  else
608  Y(upd_xpay_double)(ths->f_hat_iter, ths->alpha_iter, ths->p_hat_iter,
609  ths->mv->N_total);
610 
611  /*-----------------*/
612  Y(upd_xpay_double)(ths->r_iter, -ths->alpha_iter, ths->v_iter,
613  ths->mv->M_total);
614 
615  if(ths->flags & PRECOMPUTE_WEIGHT)
616  ths->dot_r_iter = Y(dot_w_double)(ths->r_iter,ths->w,ths->mv->M_total);
617  else
618  ths->dot_r_iter = Y(dot_double)(ths->r_iter, ths->mv->M_total);
619 
620  /*-----------------*/
621  if(ths->flags & PRECOMPUTE_WEIGHT)
622  Y(cp_w_double)(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
623  else
624  Y(cp_double)(ths->mv->f, ths->r_iter, ths->mv->M_total);
625 
626  RSWAP(ths->z_hat_iter,ths->mv->f_hat);
627  ths->mv->mv_adjoint(ths->mv);
628  RSWAP(ths->z_hat_iter,ths->mv->f_hat);
629 
630  ths->dot_z_hat_iter_old = ths->dot_z_hat_iter;
631  if(ths->flags & PRECOMPUTE_DAMP)
632  ths->dot_z_hat_iter = Y(dot_w_double)(ths->z_hat_iter, ths->w_hat,
633  ths->mv->N_total);
634  else
635  ths->dot_z_hat_iter = Y(dot_double)(ths->z_hat_iter, ths->mv->N_total);
636 
637  /*-----------------*/
638  ths->beta_iter = ths->dot_z_hat_iter / ths->dot_z_hat_iter_old;
639 
640  /*-----------------*/
641  Y(upd_axpy_double)(ths->p_hat_iter, ths->beta_iter, ths->z_hat_iter,
642  ths->mv->N_total);
643 } /* void solver_loop_one_step_cgnr */
644 
646 static void solver_loop_one_step_cgne_double(X(plan_double) *ths)
647 {
648  ths->alpha_iter = ths->dot_r_iter / ths->dot_p_hat_iter;
649 
650  /*-----------------*/
651  if(ths->flags & PRECOMPUTE_DAMP)
652  Y(upd_xpawy_double)(ths->f_hat_iter, ths->alpha_iter, ths->w_hat,
653  ths->p_hat_iter, ths->mv->N_total);
654  else
655  Y(upd_xpay_double)(ths->f_hat_iter, ths->alpha_iter, ths->p_hat_iter,
656  ths->mv->N_total);
657 
658  /*-----------------*/
659  if(ths->flags & PRECOMPUTE_DAMP)
660  Y(cp_w_double)(ths->mv->f_hat, ths->w_hat, ths->p_hat_iter,
661  ths->mv->N_total);
662  else
663  Y(cp_double)(ths->mv->f_hat, ths->p_hat_iter, ths->mv->N_total);
664 
665  ths->mv->mv_trafo(ths->mv);
666 
667  Y(upd_xpay_double)(ths->r_iter, -ths->alpha_iter, ths->mv->f,
668  ths->mv->M_total);
669 
670  ths->dot_r_iter_old = ths->dot_r_iter;
671  if(ths->flags & PRECOMPUTE_WEIGHT)
672  ths->dot_r_iter = Y(dot_w_double)(ths->r_iter,ths->w,ths->mv->M_total);
673  else
674  ths->dot_r_iter = Y(dot_double)(ths->r_iter, ths->mv->M_total);
675 
676  /*-----------------*/
677  ths->beta_iter = ths->dot_r_iter / ths->dot_r_iter_old;
678 
679  /*-----------------*/
680  if(ths->flags & PRECOMPUTE_WEIGHT)
681  Y(cp_w_double)(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total);
682  else
683  Y(cp_double)(ths->mv->f, ths->r_iter, ths->mv->M_total);
684 
685  ths->mv->mv_adjoint(ths->mv);
686 
687  Y(upd_axpy_double)(ths->p_hat_iter, ths->beta_iter, ths->mv->f_hat,
688  ths->mv->N_total);
689 
690  if(ths->flags & PRECOMPUTE_DAMP)
691  ths->dot_p_hat_iter = Y(dot_w_double)(ths->p_hat_iter, ths->w_hat,
692  ths->mv->N_total);
693  else
694  ths->dot_p_hat_iter = Y(dot_double)(ths->p_hat_iter, ths->mv->N_total);
695 } /* void solver_loop_one_step_cgne */
696 
698 void X(loop_one_step_double)(X(plan_double) *ths)
699 {
700  if(ths->flags & LANDWEBER)
701  solver_loop_one_step_landweber_double(ths);
702 
703  if(ths->flags & STEEPEST_DESCENT)
704  solver_loop_one_step_steepest_descent_double(ths);
705 
706  if(ths->flags & CGNR)
707  solver_loop_one_step_cgnr_double(ths);
708 
709  if(ths->flags & CGNE)
710  solver_loop_one_step_cgne_double(ths);
711 } /* void solver_loop_one_step */
712 
714 void X(finalize_double)(X(plan_double) *ths)
715 {
716  if(ths->flags & PRECOMPUTE_WEIGHT)
717  Y(free)(ths->w);
718 
719  if(ths->flags & PRECOMPUTE_DAMP)
720  Y(free)(ths->w_hat);
721 
722  if(ths->flags & CGNR)
723  {
724  Y(free)(ths->v_iter);
725  Y(free)(ths->z_hat_iter);
726  }
727 
728  if(ths->flags & STEEPEST_DESCENT)
729  Y(free)(ths->v_iter);
730 
731  Y(free)(ths->p_hat_iter);
732  Y(free)(ths->f_hat_iter);
733 
734  Y(free)(ths->r_iter);
735  Y(free)(ths->y);
736 }
PRECOMPUTE_DAMP
#define PRECOMPUTE_DAMP
Definition: nfft3.h:792
PRECOMPUTE_WEIGHT
#define PRECOMPUTE_WEIGHT
Definition: nfft3.h:791
STEEPEST_DESCENT
#define STEEPEST_DESCENT
Definition: nfft3.h:787
solver_loop_one_step_cgnr_double
static void solver_loop_one_step_cgnr_double(CONCAT(solver_, plan_double)*ths)
void solver_loop_one_step_cgnr
Definition: solver.c:583
CGNE
#define CGNE
Definition: nfft3.h:789
solver_loop_one_step_cgnr_complex
static void solver_loop_one_step_cgnr_complex(CONCAT(solver_, plan_complex)*ths)
void solver_loop_one_step_cgnr
Definition: solver.c:232
LANDWEBER
#define LANDWEBER
Definition: nfft3.h:786
solver_loop_one_step_cgne_double
static void solver_loop_one_step_cgne_double(CONCAT(solver_, plan_double)*ths)
void solver_loop_one_step_cgne
Definition: solver.c:646
solver_loop_one_step_landweber_complex
static void solver_loop_one_step_landweber_complex(CONCAT(solver_, plan_complex)*ths)
void solver_loop_one_step_landweber
Definition: solver.c:128
X
#define X(name)
Include header for C99 complex datatype.
Definition: fastsum.h:57
RSWAP
#define RSWAP(x, y)
Swap two vectors.
Definition: infft.h:143
solver_loop_one_step_steepest_descent_double
static void solver_loop_one_step_steepest_descent_double(CONCAT(solver_, plan_double)*ths)
void solver_loop_one_step_steepest_descent
Definition: solver.c:528
CGNR
#define CGNR
Definition: nfft3.h:788
solver_loop_one_step_landweber_double
static void solver_loop_one_step_landweber_double(CONCAT(solver_, plan_double)*ths)
void solver_loop_one_step_landweber
Definition: solver.c:479
solver_loop_one_step_cgne_complex
static void solver_loop_one_step_cgne_complex(CONCAT(solver_, plan_complex)*ths)
void solver_loop_one_step_cgne
Definition: solver.c:295
NORMS_FOR_LANDWEBER
#define NORMS_FOR_LANDWEBER
Definition: nfft3.h:790
solver_loop_one_step_steepest_descent_complex
static void solver_loop_one_step_steepest_descent_complex(CONCAT(solver_, plan_complex)*ths)
void solver_loop_one_step_steepest_descent
Definition: solver.c:177
nfft3.h
Header file for the nfft3 library.
CSWAP
#define CSWAP(x, y)
Swap two vectors.
Definition: infft.h:139
Generated on Thu Apr 12 2018 21:00:50 for NFFT by   doxygen 1.8.6