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