NFFT  3.3.0
glacier.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 
39 static NFFT_R my_weight(NFFT_R z, NFFT_R a, NFFT_R b, NFFT_R c)
40 {
41  return NFFT_M(pow)(NFFT_K(0.25) - z * z, b) / (c + NFFT_M(pow)(NFFT_M(fabs)(z), NFFT_K(2.0) * a));
42 }
43 
45 static void glacier(int N, int M)
46 {
47  int j, k, k0, k1, l, my_N[2], my_n[2];
48  NFFT_R tmp_y;
49  NFFT(plan) p;
50  SOLVER(plan_complex) ip;
51  FILE* fp;
52 
53  /* initialise p */
54  my_N[0] = N;
55  my_n[0] = (int)NFFT(next_power_of_2)(N);
56  my_N[1] = N;
57  my_n[1] = (int)NFFT(next_power_of_2)(N);
58  NFFT(init_guru)(&p, 2, my_N, M, my_n, 6,
59  PRE_PHI_HUT | PRE_FULL_PSI |
60  MALLOC_X | MALLOC_F_HAT | MALLOC_F |
61  FFTW_INIT | FFT_OUT_OF_PLACE,
62  FFTW_MEASURE | FFTW_DESTROY_INPUT);
63 
64  /* initialise ip, specific */
65  SOLVER(init_advanced_complex)(&ip, (NFFT(mv_plan_complex)*) (&p),
66  CGNE | PRECOMPUTE_DAMP);
67  fprintf(stderr, "Using the generic solver!");
68 
69  /* init nodes */
70  fp = fopen("input_data.dat", "r");
71  for (j = 0; j < p.M_total; j++)
72  {
73  fscanf(fp, "%" NFFT__FES__ " %" NFFT__FES__ " %" NFFT__FES__, &p.x[2 * j + 0], &p.x[2 * j + 1], &tmp_y);
74  ip.y[j] = tmp_y;
75  }
76  fclose(fp);
77 
78  /* precompute psi */
79  if (p.flags & PRE_ONE_PSI)
80  NFFT(precompute_one_psi)(&p);
81 
82  /* initialise damping factors */
83  if (ip.flags & PRECOMPUTE_DAMP)
84  for (k0 = 0; k0 < p.N[0]; k0++)
85  for (k1 = 0; k1 < p.N[1]; k1++)
86  ip.w_hat[k0 * p.N[1] + k1] = my_weight(((NFFT_R) ((NFFT_R)(k0) - (NFFT_R)(p.N[0]) / NFFT_K(2.0))) / ((NFFT_R)(p.N[0])),
87  NFFT_K(0.5), NFFT_K(3.0), NFFT_K(0.001))
88  * my_weight((((NFFT_R)(k1) - (NFFT_R)(p.N[1]) / NFFT_K(2.0))) / ((NFFT_R)(p.N[1])), NFFT_K(0.5), NFFT_K(3.0), NFFT_K(0.001));
89 
90  /* init some guess */
91  for (k = 0; k < p.N_total; k++)
92  ip.f_hat_iter[k] = NFFT_K(0.0);
93 
94  /* inverse trafo */
95  SOLVER(before_loop_complex)(&ip);
96 
97  for (l = 0; l < 40; l++)
98  {
99  fprintf(stderr, "Residual ||r||=%" NFFT__FES__ ",\n", NFFT_M(sqrt)(ip.dot_r_iter));
100  SOLVER(loop_one_step_complex)(&ip);
101  }
102 
103  for (k = 0; k < p.N_total; k++)
104  printf("%" NFFT__FES__ " %" NFFT__FES__ "\n", NFFT_M(creal)(ip.f_hat_iter[k]), NFFT_M(cimag)(ip.f_hat_iter[k]));
105 
106  SOLVER(finalize_complex)(&ip);
107  NFFT(finalize)(&p);
108 }
109 
111 static void glacier_cv(int N, int M, int M_cv, unsigned solver_flags)
112 {
113  int j, k, k0, k1, l, my_N[2], my_n[2];
114  NFFT_R tmp_y, r;
115  NFFT(plan) p, cp;
116  SOLVER(plan_complex) ip;
117  NFFT_C* cp_y;
118  FILE* fp;
119  int M_re = M - M_cv;
120 
121  /* initialise p for reconstruction */
122  my_N[0] = N;
123  my_n[0] = (int)NFFT(next_power_of_2)(N);
124  my_N[1] = N;
125  my_n[1] = (int)NFFT(next_power_of_2)(N);
126  NFFT(init_guru)(&p, 2, my_N, M_re, my_n, 6,
127  PRE_PHI_HUT | PRE_FULL_PSI |
128  MALLOC_X | MALLOC_F_HAT | MALLOC_F |
129  FFTW_INIT | FFT_OUT_OF_PLACE,
130  FFTW_MEASURE | FFTW_DESTROY_INPUT);
131 
132  /* initialise ip, specific */
133  SOLVER(init_advanced_complex)(&ip, (NFFT(mv_plan_complex)*) (&p), solver_flags);
134 
135  /* initialise cp for validation */
136  cp_y = (NFFT_C*) NFFT(malloc)((size_t)(M) * sizeof(NFFT_C));
137  NFFT(init_guru)(&cp, 2, my_N, M, my_n, 6,
138  PRE_PHI_HUT | PRE_FULL_PSI |
139  MALLOC_X | MALLOC_F |
140  FFTW_INIT | FFT_OUT_OF_PLACE,
141  FFTW_MEASURE | FFTW_DESTROY_INPUT);
142 
143  cp.f_hat = ip.f_hat_iter;
144 
145  /* set up data in cp and cp_y */
146  fp = fopen("input_data.dat", "r");
147  for (j = 0; j < cp.M_total; j++)
148  {
149  fscanf(fp, "%" NFFT__FES__ " %" NFFT__FES__ " %" NFFT__FES__, &cp.x[2 * j + 0], &cp.x[2 * j + 1], &tmp_y);
150  cp_y[j] = tmp_y;
151  }
152  fclose(fp);
153 
154  /* copy part of the data to p and ip */
155  for (j = 0; j < p.M_total; j++)
156  {
157  p.x[2 * j + 0] = cp.x[2 * j + 0];
158  p.x[2 * j + 1] = cp.x[2 * j + 1];
159  ip.y[j] = tmp_y;
160  }
161 
162  /* precompute psi */
163  if (p.flags & PRE_ONE_PSI)
164  NFFT(precompute_one_psi)(&p);
165 
166  /* precompute psi */
167  if (cp.flags & PRE_ONE_PSI)
168  NFFT(precompute_one_psi)(&cp);
169 
170  /* initialise damping factors */
171  if (ip.flags & PRECOMPUTE_DAMP)
172  for (k0 = 0; k0 < p.N[0]; k0++)
173  for (k1 = 0; k1 < p.N[1]; k1++)
174  ip.w_hat[k0 * p.N[1] + k1] = my_weight((((NFFT_R)(k0) - (NFFT_R)(p.N[0]) / NFFT_K(2.0))) / ((NFFT_R)(p.N[0])),
175  NFFT_K(0.5), NFFT_K(3.0), NFFT_K(0.001))
176  * my_weight((((NFFT_R)(k1) - (NFFT_R)(p.N[1]) / NFFT_K(2.0))) / ((NFFT_R)(p.N[1])), NFFT_K(0.5), NFFT_K(3.0), NFFT_K(0.001));
177 
178  /* init some guess */
179  for (k = 0; k < p.N_total; k++)
180  ip.f_hat_iter[k] = NFFT_K(0.0);
181 
182  /* inverse trafo */
183  SOLVER(before_loop_complex)(&ip);
184  // fprintf(stderr,"iteration starts,\t");
185  for (l = 0; l < 40; l++)
186  SOLVER(loop_one_step_complex)(&ip);
187 
188  //fprintf(stderr,"r=%1.2e, ",sqrt(ip.dot_r_iter)/M_re);
189 
190  NFFT_CSWAP(p.f_hat, ip.f_hat_iter);
191  NFFT(trafo)(&p);
192  NFFT_CSWAP(p.f_hat, ip.f_hat_iter);
193  NFFT(upd_axpy_complex)(p.f, -1, ip.y, M_re);
194  r = NFFT_M(sqrt)(NFFT(dot_complex)(p.f, M_re) / NFFT(dot_complex)(cp_y, M));
195  fprintf(stderr, "r=%1.2" NFFT__FES__ ", ", r);
196  printf("$%1.1" NFFT__FES__ "$ & ", r);
197 
198  NFFT(trafo)(&cp);
199  NFFT(upd_axpy_complex)(&cp.f[M_re], -1, &cp_y[M_re], M_cv);
200  r = NFFT_M(sqrt)(NFFT(dot_complex)(&cp.f[M_re], M_cv) / NFFT(dot_complex)(cp_y, M));
201  fprintf(stderr, "r_1=%1.2" NFFT__FES__ "\t", r);
202  printf("$%1.1" NFFT__FES__ "$ & ", r);
203 
204  NFFT(finalize)(&cp);
205  SOLVER(finalize_complex)(&ip);
206  NFFT(finalize)(&p);
207 }
208 
210 int main(int argc, char **argv)
211 {
212  int M_cv;
213 
214  if (argc < 3)
215  {
216  fprintf(stderr, "Call this program from the Matlab script glacier.m!");
217  return EXIT_FAILURE;
218  }
219 
220  if (argc == 3)
221  glacier(atoi(argv[1]), atoi(argv[2]));
222  else
223  for (M_cv = atoi(argv[3]); M_cv <= atoi(argv[5]); M_cv += atoi(argv[4]))
224  {
225  fprintf(stderr, "\nM_cv=%d,\t", M_cv);
226  printf("$%d$ & ", M_cv);
227  fprintf(stderr, "cgne+damp: ");
228  glacier_cv(atoi(argv[1]), atoi(argv[2]), M_cv, CGNE | PRECOMPUTE_DAMP);
229  //fprintf(stderr,"cgne: ");
230  //glacier_cv(atoi(argv[1]),atoi(argv[2]),M_cv,CGNE);
231  fprintf(stderr, "cgnr: ");
232  glacier_cv(atoi(argv[1]), atoi(argv[2]), M_cv, CGNR);
233  fprintf(stderr, "cgnr: ");
234  glacier_cv(atoi(argv[1]) / 4, atoi(argv[2]), M_cv, CGNR);
235  printf("XXX \\\\\n");
236  }
237 
238  fprintf(stderr, "\n");
239 
240  return EXIT_SUCCESS;
241 }
242 /* \} */
static NFFT_R my_weight(NFFT_R z, NFFT_R a, NFFT_R b, NFFT_R c)
Generalised Sobolev weight.
Definition: glacier.c:39
static void glacier(int N, int M)
Reconstruction routine.
Definition: glacier.c:45
static void glacier_cv(int N, int M, int M_cv, unsigned solver_flags)
Reconstruction routine with cross validation.
Definition: glacier.c:111