NFFT 3.5.3alpha
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
40void 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
76void X(init_complex)(X(plan_complex)* ths, Y(mv_plan_complex) *mv)
77{
78 X(init_advanced_complex)(ths, mv, CGNR);
79}
80
81void 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
128static 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
177static 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
232static 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
295static 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
347void X(loop_one_step_complex)(X(plan_complex) *ths)
348{
349 if(ths->flags & LANDWEBER)
351
352 if(ths->flags & STEEPEST_DESCENT)
354
355 if(ths->flags & CGNR)
357
358 if(ths->flags & CGNE)
360} /* void solver_loop_one_step */
361
363void 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
392void 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
427void X(init_double)(X(plan_double)* ths, Y(mv_plan_double) *mv)
428{
429 X(init_advanced_double)(ths, mv, CGNR);
430}
431
432void 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
479static 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
528static 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
583static 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
646static 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
698void X(loop_one_step_double)(X(plan_double) *ths)
699{
700 if(ths->flags & LANDWEBER)
702
703 if(ths->flags & STEEPEST_DESCENT)
705
706 if(ths->flags & CGNR)
708
709 if(ths->flags & CGNE)
711} /* void solver_loop_one_step */
712
714void 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}
#define RSWAP(x, y)
Swap two vectors.
Definition infft.h:143
#define CSWAP(x, y)
Swap two vectors.
Definition infft.h:139
#define STEEPEST_DESCENT
Definition nfft3.h:807
#define CGNR
Definition nfft3.h:808
#define NORMS_FOR_LANDWEBER
Definition nfft3.h:810
#define LANDWEBER
Definition nfft3.h:806
#define PRECOMPUTE_DAMP
Definition nfft3.h:812
#define PRECOMPUTE_WEIGHT
Definition nfft3.h:811
#define CGNE
Definition nfft3.h:809
Internal header file for auxiliary definitions and functions.
Header file for the nfft3 library.
static void solver_loop_one_step_landweber_complex(CONCAT(solver_, plan_complex) *ths)
void solver_loop_one_step_landweber
Definition solver.c:128
static void solver_loop_one_step_cgne_complex(CONCAT(solver_, plan_complex) *ths)
void solver_loop_one_step_cgne
Definition solver.c:295
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
static void solver_loop_one_step_cgnr_double(CONCAT(solver_, plan_double) *ths)
void solver_loop_one_step_cgnr
Definition solver.c:583
static void solver_loop_one_step_landweber_double(CONCAT(solver_, plan_double) *ths)
void solver_loop_one_step_landweber
Definition solver.c:479
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
static void solver_loop_one_step_cgnr_complex(CONCAT(solver_, plan_complex) *ths)
void solver_loop_one_step_cgnr
Definition solver.c:232
static void solver_loop_one_step_cgne_double(CONCAT(solver_, plan_double) *ths)
void solver_loop_one_step_cgne
Definition solver.c:646