NFFT  3.3.0
solver/simple_test.c
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 
21 #include <stdio.h>
22 #include <math.h>
23 #include <string.h>
24 #include <stdlib.h>
25 #include <complex.h>
26 
27 #define NFFT_PRECISION_DOUBLE
28 
29 #include "nfft3mp.h"
30 
31 /* void simple_test_infft_1d(int N, int M, int iter) */
32 /* { */
33 /* int k,l; /\**< index for nodes, freqencies,iter*\/ */
34 /* nfft_plan p; /\**< plan for the nfft *\/ */
35 /* infft_plan ip; /\**< plan for the inverse nfft *\/ */
36 
37 /* /\** initialise an one dimensional plan *\/ */
38 /* nfft_init_1d(&p, N, M); */
39 
40 /* /\** init pseudo random nodes *\/ */
41 /* nfft_vrand_shifted_unit_double(p.x,p.M_total); */
42 
43 /* /\** precompute psi, the entries of the matrix B *\/ */
44 /* if(p.flags & PRE_ONE_PSI) */
45 /* nfft_precompute_one_psi(&p); */
46 
47 /* /\** initialise inverse plan *\/ */
48 /* infft_init(&ip,&p); */
49 
50 /* /\** init pseudo random samples and show them *\/ */
51 /* nfft_vrand_unit_complex(ip.y,p.M_total); */
52 /* nfft_vpr_complex(ip.y,p.M_total,"Given data, vector y"); */
53 
54 /* /\** initialise some guess f_hat_0 and solve *\/ */
55 /* for(k=0;k<p.N_total;k++) */
56 /* ip.f_hat_iter[k]=0; */
57 
58 /* nfft_vpr_complex(ip.f_hat_iter,p.N_total,"Initial guess, vector f_hat_iter"); */
59 
60 /* CSWAP(ip.f_hat_iter,p.f_hat); */
61 /* nfft_trafo(&p); */
62 /* nfft_vpr_complex(p.f,p.M_total,"Data fit, vector f"); */
63 /* CSWAP(ip.f_hat_iter,p.f_hat); */
64 
65 /* infft_before_loop(&ip); */
66 /* printf("\n Residual r=%e\n",ip.dot_r_iter); */
67 
68 /* for(l=0;l<iter;l++) */
69 /* { */
70 /* printf("\n********** Iteration l=%d **********\n",l); */
71 /* infft_loop_one_step(&ip); */
72 /* nfft_vpr_complex(ip.f_hat_iter,p.N_total, */
73 /* "Approximate solution, vector f_hat_iter"); */
74 
75 /* CSWAP(ip.f_hat_iter,p.f_hat); */
76 /* nfft_trafo(&p); */
77 /* nfft_vpr_complex(p.f,p.M_total,"Data fit, vector f"); */
78 /* CSWAP(ip.f_hat_iter,p.f_hat); */
79 
80 /* printf("\n Residual r=%e\n",ip.dot_r_iter); */
81 /* } */
82 
83 /* infft_finalize(&ip); */
84 /* nfft_finalize(&p); */
85 /* } */
86 
88 static void simple_test_solver_nfft_1d(int N, int M, int iter)
89 {
90  int k, l;
91  NFFT(plan) p;
92  SOLVER(plan_complex) ip;
93  const char *error_str;
94 
96  NFFT(init_1d)(&p, N, M);
97 
99  NFFT(vrand_shifted_unit_double)(p.x, p.M_total);
100 
102  if (p.flags & PRE_ONE_PSI)
103  NFFT(precompute_one_psi)(&p);
104 
106  SOLVER(init_complex)(&ip, (NFFT(mv_plan_complex)*) (&p));
107 
109  NFFT(vrand_unit_complex)(ip.y, p.M_total);
110  NFFT(vpr_complex)(ip.y, p.M_total, "Given data, vector y");
111 
113  for (k = 0; k < p.N_total; k++)
114  ip.f_hat_iter[k] = NFFT_K(0.0);
115 
116  NFFT(vpr_complex)(ip.f_hat_iter, p.N_total,
117  "Initial guess, vector f_hat_iter");
118 
120  error_str = NFFT(check)(&p);
121  if (error_str != 0)
122  {
123  printf("Error in nfft module: %s\n", error_str);
124  return;
125  }
126 
127  NFFT_CSWAP(ip.f_hat_iter, p.f_hat);
128  NFFT(trafo)(&p);
129  NFFT(vpr_complex)(p.f, p.M_total, "Data fit, vector f");
130  NFFT_CSWAP(ip.f_hat_iter, p.f_hat);
131 
132  SOLVER(before_loop_complex)(&ip);
133  printf("\n Residual r=%" NFFT__FES__ "\n", ip.dot_r_iter);
134 
135  for (l = 0; l < iter; l++)
136  {
137  printf("\n********** Iteration l=%d **********\n", l);
138  SOLVER(loop_one_step_complex)(&ip);
139  NFFT(vpr_complex)(ip.f_hat_iter, p.N_total,
140  "Approximate solution, vector f_hat_iter");
141 
142  NFFT_CSWAP(ip.f_hat_iter, p.f_hat);
143  NFFT(trafo)(&p);
144  NFFT(vpr_complex)(p.f, p.M_total, "Data fit, vector f");
145  NFFT_CSWAP(ip.f_hat_iter, p.f_hat);
146 
147  printf("\n Residual r=%" NFFT__FES__ "\n", ip.dot_r_iter);
148  }
149 
150  SOLVER(finalize_complex)(&ip);
151  NFFT(finalize)(&p);
152 }
153 
155 int main(void)
156 {
157  printf("\n Computing a one dimensional inverse nfft\n");
158 
159  simple_test_solver_nfft_1d(16, 8, 9);
160 
161  return EXIT_SUCCESS;
162 }
163 /* \} */