00001
00007 #include "util.h"
00008 #include "nfft3.h"
00009
00021 #define F(MV, FLT, FLT_TYPE, NAME, ...) void i ## MV ## _ ## NAME(__VA_ARGS__)
00022
00023
00024 #define MACRO_SOLVER_IMPL(MV, FLT, FLT_TYPE) \
00025 \
00026 F(MV, FLT, FLT_TYPE, init_advanced, i ## MV ## _plan *ths, MV ## _plan *mv, \
00027 unsigned i ## MV ## _flags) \
00028 { \
00029 ths->mv = mv; \
00030 ths->flags = i ## MV ## _flags; \
00031 \
00032 ths->y = (FLT_TYPE*)fftw_malloc(ths->mv->M_total*sizeof(FLT_TYPE));\
00033 ths->r_iter = (FLT_TYPE*)fftw_malloc(ths->mv->M_total*sizeof(FLT_TYPE));\
00034 ths->f_hat_iter = (FLT_TYPE*)fftw_malloc(ths->mv->N_total*sizeof(FLT_TYPE));\
00035 ths->p_hat_iter = (FLT_TYPE*)fftw_malloc(ths->mv->N_total*sizeof(FLT_TYPE));\
00036 \
00037 if(ths->flags & LANDWEBER) \
00038 ths->z_hat_iter = ths->p_hat_iter; \
00039 \
00040 if(ths->flags & STEEPEST_DESCENT) \
00041 { \
00042 ths->z_hat_iter = ths->p_hat_iter; \
00043 ths->v_iter = (FLT_TYPE*)fftw_malloc(ths->mv->M_total* \
00044 sizeof(FLT_TYPE)); \
00045 } \
00046 \
00047 if(ths->flags & CGNR) \
00048 { \
00049 ths->z_hat_iter = (FLT_TYPE*)fftw_malloc(ths->mv->N_total* \
00050 sizeof(FLT_TYPE)); \
00051 ths->v_iter = (FLT_TYPE*)fftw_malloc(ths->mv->M_total* \
00052 sizeof(FLT_TYPE)); \
00053 } \
00054 \
00055 if(ths->flags & CGNE) \
00056 ths->z_hat_iter = ths->p_hat_iter; \
00057 \
00058 if(ths->flags & PRECOMPUTE_WEIGHT) \
00059 ths->w = (double*) fftw_malloc(ths->mv->M_total*sizeof(double)); \
00060 \
00061 if(ths->flags & PRECOMPUTE_DAMP) \
00062 ths->w_hat = (double*) fftw_malloc(ths->mv->N_total*sizeof(double)); \
00063 } \
00064 \
00065 \
00066 F(MV, FLT, FLT_TYPE, init, i ## MV ## _plan *ths, MV ## _plan *mv) \
00067 { \
00068 i ## MV ## _init_advanced(ths, mv, CGNR); \
00069 } \
00070 \
00071 \
00072 F(MV, FLT, FLT_TYPE, before_loop, i ## MV ## _plan *ths) \
00073 { \
00074 nfft_cp_ ## FLT(ths->mv->f_hat, ths->f_hat_iter, ths->mv->N_total); \
00075 \
00076 NFFT_SWAP_ ## FLT(ths->r_iter, ths->mv->f); \
00077 MV ## _trafo(ths->mv); \
00078 NFFT_SWAP_ ## FLT(ths->r_iter, ths->mv->f); \
00079 \
00080 nfft_upd_axpy_ ## FLT(ths->r_iter, -1.0, ths->y, ths->mv->M_total); \
00081 \
00082 if((!(ths->flags & LANDWEBER)) || (ths->flags & NORMS_FOR_LANDWEBER)) \
00083 { \
00084 if(ths->flags & PRECOMPUTE_WEIGHT) \
00085 ths->dot_r_iter = \
00086 nfft_dot_w_ ## FLT(ths->r_iter, ths->w, ths->mv->M_total); \
00087 else \
00088 ths->dot_r_iter = nfft_dot_ ## FLT(ths->r_iter, ths->mv->M_total); \
00089 } \
00090 \
00091 \
00092 if(ths->flags & PRECOMPUTE_WEIGHT) \
00093 nfft_cp_w_ ## FLT(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total); \
00094 else \
00095 nfft_cp_ ## FLT(ths->mv->f, ths->r_iter, ths->mv->M_total); \
00096 \
00097 NFFT_SWAP_ ## FLT(ths->z_hat_iter, ths->mv->f_hat); \
00098 MV ## _adjoint(ths->mv); \
00099 NFFT_SWAP_ ## FLT(ths->z_hat_iter, ths->mv->f_hat); \
00100 \
00101 if((!(ths->flags & LANDWEBER)) || (ths->flags & NORMS_FOR_LANDWEBER)) \
00102 { \
00103 if(ths->flags & PRECOMPUTE_DAMP) \
00104 ths->dot_z_hat_iter = \
00105 nfft_dot_w_ ## FLT(ths->z_hat_iter, ths->w_hat, ths->mv->N_total); \
00106 else \
00107 ths->dot_z_hat_iter = nfft_dot_ ## FLT(ths->z_hat_iter, ths->mv->N_total); \
00108 } \
00109 \
00110 if(ths->flags & CGNE) \
00111 ths->dot_p_hat_iter = ths->dot_z_hat_iter; \
00112 \
00113 if(ths->flags & CGNR) \
00114 nfft_cp_ ## FLT(ths->p_hat_iter, ths->z_hat_iter, ths->mv->N_total); \
00115 } \
00116 \
00117 \
00118 F(MV, FLT, FLT_TYPE, loop_one_step_landweber, i ## MV ## _plan *ths) \
00119 { \
00120 if(ths->flags & PRECOMPUTE_DAMP) \
00121 nfft_upd_xpawy_ ## FLT(ths->f_hat_iter, ths->alpha_iter, ths->w_hat, \
00122 ths->z_hat_iter, ths->mv->N_total); \
00123 else \
00124 nfft_upd_xpay_ ## FLT(ths->f_hat_iter, ths->alpha_iter, ths->z_hat_iter, \
00125 ths->mv->N_total); \
00126 \
00127 \
00128 nfft_cp_ ## FLT(ths->mv->f_hat, ths->f_hat_iter, ths->mv->N_total); \
00129 \
00130 NFFT_SWAP_ ## FLT(ths->r_iter,ths->mv->f); \
00131 MV ## _trafo(ths->mv); \
00132 NFFT_SWAP_ ## FLT(ths->r_iter,ths->mv->f); \
00133 \
00134 nfft_upd_axpy_ ## FLT(ths->r_iter, -1.0, ths->y, ths->mv->M_total); \
00135 \
00136 if(ths->flags & NORMS_FOR_LANDWEBER) \
00137 { \
00138 if(ths->flags & PRECOMPUTE_WEIGHT) \
00139 ths->dot_r_iter = nfft_dot_w_ ## FLT(ths->r_iter,ths->w, ths->mv->M_total); \
00140 else \
00141 ths->dot_r_iter = nfft_dot_ ## FLT(ths->r_iter, ths->mv->M_total); \
00142 } \
00143 \
00144 \
00145 if(ths->flags & PRECOMPUTE_WEIGHT) \
00146 nfft_cp_w_ ## FLT(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total); \
00147 else \
00148 nfft_cp_ ## FLT(ths->mv->f, ths->r_iter, ths->mv->M_total); \
00149 \
00150 NFFT_SWAP_ ## FLT(ths->z_hat_iter,ths->mv->f_hat); \
00151 MV ## _adjoint(ths->mv); \
00152 NFFT_SWAP_ ## FLT(ths->z_hat_iter,ths->mv->f_hat); \
00153 \
00154 if(ths->flags & NORMS_FOR_LANDWEBER) \
00155 { \
00156 if(ths->flags & PRECOMPUTE_DAMP) \
00157 ths->dot_z_hat_iter = \
00158 nfft_dot_w_ ## FLT(ths->z_hat_iter, ths->w_hat, ths->mv->N_total); \
00159 else \
00160 ths->dot_z_hat_iter = nfft_dot_ ## FLT(ths->z_hat_iter, ths->mv->N_total); \
00161 } \
00162 } \
00163 \
00164 \
00165 F(MV, FLT, FLT_TYPE, loop_one_step_steepest_descent, i ## MV ## _plan *ths) \
00166 { \
00167 if(ths->flags & PRECOMPUTE_DAMP) \
00168 nfft_cp_w_ ## FLT(ths->mv->f_hat, ths->w_hat, ths->z_hat_iter, \
00169 ths->mv->N_total); \
00170 else \
00171 nfft_cp_ ## FLT(ths->mv->f_hat, ths->z_hat_iter, ths->mv->N_total); \
00172 \
00173 NFFT_SWAP_ ## FLT(ths->v_iter,ths->mv->f); \
00174 MV ## _trafo(ths->mv); \
00175 NFFT_SWAP_ ## FLT(ths->v_iter,ths->mv->f); \
00176 \
00177 if(ths->flags & PRECOMPUTE_WEIGHT) \
00178 ths->dot_v_iter = nfft_dot_w_ ## FLT(ths->v_iter, ths->w, \
00179 ths->mv->M_total); \
00180 else \
00181 ths->dot_v_iter = nfft_dot_ ## FLT(ths->v_iter, ths->mv->M_total); \
00182 \
00183 \
00184 ths->alpha_iter = ths->dot_z_hat_iter / ths->dot_v_iter; \
00185 \
00186 \
00187 if(ths->flags & PRECOMPUTE_DAMP) \
00188 nfft_upd_xpawy_ ## FLT(ths->f_hat_iter, ths->alpha_iter, ths->w_hat, \
00189 ths->z_hat_iter, ths->mv->N_total); \
00190 else \
00191 nfft_upd_xpay_ ## FLT(ths->f_hat_iter, ths->alpha_iter, ths->z_hat_iter, \
00192 ths->mv->N_total); \
00193 \
00194 \
00195 nfft_upd_xpay_ ## FLT(ths->r_iter, -ths->alpha_iter, ths->v_iter, \
00196 ths->mv->M_total); \
00197 \
00198 if(ths->flags & PRECOMPUTE_WEIGHT) \
00199 ths->dot_r_iter = nfft_dot_w_ ## FLT(ths->r_iter, ths->w, \
00200 ths->mv->M_total); \
00201 else \
00202 ths->dot_r_iter = nfft_dot_ ## FLT(ths->r_iter, ths->mv->M_total); \
00203 \
00204 \
00205 if(ths->flags & PRECOMPUTE_WEIGHT) \
00206 nfft_cp_w_ ## FLT(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total); \
00207 else \
00208 nfft_cp_ ## FLT(ths->mv->f, ths->r_iter, ths->mv->M_total); \
00209 \
00210 NFFT_SWAP_ ## FLT(ths->z_hat_iter,ths->mv->f_hat); \
00211 MV ## _adjoint(ths->mv); \
00212 NFFT_SWAP_ ## FLT(ths->z_hat_iter,ths->mv->f_hat); \
00213 \
00214 if(ths->flags & PRECOMPUTE_DAMP) \
00215 ths->dot_z_hat_iter = \
00216 nfft_dot_w_ ## FLT(ths->z_hat_iter, ths->w_hat, ths->mv->N_total); \
00217 else \
00218 ths->dot_z_hat_iter = nfft_dot_ ## FLT(ths->z_hat_iter, ths->mv->N_total);\
00219 } \
00220 \
00221 \
00222 F(MV, FLT, FLT_TYPE, loop_one_step_cgnr, i ## MV ## _plan *ths) \
00223 { \
00224 if(ths->flags & PRECOMPUTE_DAMP) \
00225 nfft_cp_w_ ## FLT(ths->mv->f_hat, ths->w_hat, ths->p_hat_iter, \
00226 ths->mv->N_total); \
00227 else \
00228 nfft_cp_ ## FLT(ths->mv->f_hat, ths->p_hat_iter, ths->mv->N_total); \
00229 \
00230 NFFT_SWAP_ ## FLT(ths->v_iter,ths->mv->f); \
00231 MV ## _trafo(ths->mv); \
00232 NFFT_SWAP_ ## FLT(ths->v_iter,ths->mv->f); \
00233 \
00234 if(ths->flags & PRECOMPUTE_WEIGHT) \
00235 ths->dot_v_iter = nfft_dot_w_ ## FLT(ths->v_iter, ths->w, \
00236 ths->mv->M_total); \
00237 else \
00238 ths->dot_v_iter = nfft_dot_ ## FLT(ths->v_iter, ths->mv->M_total); \
00239 \
00240 \
00241 ths->alpha_iter = ths->dot_z_hat_iter / ths->dot_v_iter; \
00242 \
00243 \
00244 if(ths->flags & PRECOMPUTE_DAMP) \
00245 nfft_upd_xpawy_ ## FLT(ths->f_hat_iter, ths->alpha_iter, ths->w_hat, \
00246 ths->p_hat_iter, ths->mv->N_total); \
00247 else \
00248 nfft_upd_xpay_ ## FLT(ths->f_hat_iter, ths->alpha_iter, ths->p_hat_iter, \
00249 ths->mv->N_total); \
00250 \
00251 \
00252 nfft_upd_xpay_ ## FLT(ths->r_iter, -ths->alpha_iter, ths->v_iter, \
00253 ths->mv->M_total); \
00254 \
00255 if(ths->flags & PRECOMPUTE_WEIGHT) \
00256 ths->dot_r_iter = nfft_dot_w_ ## FLT(ths->r_iter, ths->w, \
00257 ths->mv->M_total); \
00258 else \
00259 ths->dot_r_iter = nfft_dot_ ## FLT(ths->r_iter, ths->mv->M_total); \
00260 \
00261 \
00262 if(ths->flags & PRECOMPUTE_WEIGHT) \
00263 nfft_cp_w_ ## FLT(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total); \
00264 else \
00265 nfft_cp_ ## FLT(ths->mv->f, ths->r_iter, ths->mv->M_total); \
00266 \
00267 NFFT_SWAP_ ## FLT(ths->z_hat_iter,ths->mv->f_hat); \
00268 MV ## _adjoint(ths->mv); \
00269 NFFT_SWAP_ ## FLT(ths->z_hat_iter,ths->mv->f_hat); \
00270 \
00271 ths->dot_z_hat_iter_old = ths->dot_z_hat_iter; \
00272 if(ths->flags & PRECOMPUTE_DAMP) \
00273 ths->dot_z_hat_iter = \
00274 nfft_dot_w_ ## FLT(ths->z_hat_iter, ths->w_hat, ths->mv->N_total); \
00275 else \
00276 ths->dot_z_hat_iter = nfft_dot_ ## FLT(ths->z_hat_iter, ths->mv->N_total);\
00277 \
00278 \
00279 ths->beta_iter = ths->dot_z_hat_iter / ths->dot_z_hat_iter_old; \
00280 \
00281 \
00282 nfft_upd_axpy_ ## FLT(ths->p_hat_iter, ths->beta_iter, ths->z_hat_iter, \
00283 ths->mv->N_total); \
00284 } \
00285 \
00286 \
00287 F(MV, FLT, FLT_TYPE, loop_one_step_cgne, i ## MV ## _plan *ths) \
00288 { \
00289 ths->alpha_iter = ths->dot_r_iter / ths->dot_p_hat_iter; \
00290 \
00291 \
00292 if(ths->flags & PRECOMPUTE_DAMP) \
00293 nfft_upd_xpawy_ ## FLT(ths->f_hat_iter, ths->alpha_iter, ths->w_hat, \
00294 ths->p_hat_iter, ths->mv->N_total); \
00295 else \
00296 nfft_upd_xpay_ ## FLT(ths->f_hat_iter, ths->alpha_iter, ths->p_hat_iter, \
00297 ths->mv->N_total); \
00298 \
00299 \
00300 if(ths->flags & PRECOMPUTE_DAMP) \
00301 nfft_cp_w_ ## FLT(ths->mv->f_hat, ths->w_hat, ths->p_hat_iter, \
00302 ths->mv->N_total); \
00303 else \
00304 nfft_cp_ ## FLT(ths->mv->f_hat, ths->p_hat_iter, ths->mv->N_total); \
00305 \
00306 MV ## _trafo(ths->mv); \
00307 \
00308 nfft_upd_xpay_ ## FLT(ths->r_iter, -ths->alpha_iter, ths->mv->f, \
00309 ths->mv->M_total); \
00310 \
00311 ths->dot_r_iter_old = ths->dot_r_iter; \
00312 if(ths->flags & PRECOMPUTE_WEIGHT) \
00313 ths->dot_r_iter = nfft_dot_w_ ## FLT(ths->r_iter, ths->w, \
00314 ths->mv->M_total); \
00315 else \
00316 ths->dot_r_iter = nfft_dot_ ## FLT(ths->r_iter, ths->mv->M_total); \
00317 \
00318 \
00319 ths->beta_iter = ths->dot_r_iter / ths->dot_r_iter_old; \
00320 \
00321 \
00322 if(ths->flags & PRECOMPUTE_WEIGHT) \
00323 nfft_cp_w_ ## FLT(ths->mv->f, ths->w, ths->r_iter, ths->mv->M_total); \
00324 else \
00325 nfft_cp_ ## FLT(ths->mv->f, ths->r_iter, ths->mv->M_total); \
00326 \
00327 MV ## _adjoint(ths->mv); \
00328 \
00329 nfft_upd_axpy_ ## FLT(ths->p_hat_iter, ths->beta_iter, ths->mv->f_hat, \
00330 ths->mv->N_total); \
00331 \
00332 if(ths->flags & PRECOMPUTE_DAMP) \
00333 ths->dot_p_hat_iter = \
00334 nfft_dot_w_ ## FLT(ths->p_hat_iter, ths->w_hat, ths->mv->N_total); \
00335 else \
00336 ths->dot_p_hat_iter = nfft_dot_ ## FLT(ths->p_hat_iter, ths->mv->N_total);\
00337 } \
00338 \
00339 \
00340 F(MV, FLT, FLT_TYPE, loop_one_step, i ## MV ## _plan *ths) \
00341 { \
00342 if(ths->flags & LANDWEBER) \
00343 i ## MV ## _loop_one_step_landweber(ths); \
00344 \
00345 if(ths->flags & STEEPEST_DESCENT) \
00346 i ## MV ## _loop_one_step_steepest_descent(ths); \
00347 \
00348 if(ths->flags & CGNR) \
00349 i ## MV ## _loop_one_step_cgnr(ths); \
00350 \
00351 if(ths->flags & CGNE) \
00352 i ## MV ## _loop_one_step_cgne(ths); \
00353 } \
00354 \
00355 \
00356 F(MV, FLT, FLT_TYPE, finalize, i ## MV ## _plan *ths) \
00357 { \
00358 if(ths->flags & PRECOMPUTE_WEIGHT) \
00359 fftw_free(ths->w); \
00360 \
00361 if(ths->flags & PRECOMPUTE_DAMP) \
00362 fftw_free(ths->w_hat); \
00363 \
00364 if(ths->flags & CGNR) \
00365 { \
00366 fftw_free(ths->v_iter); \
00367 fftw_free(ths->z_hat_iter); \
00368 } \
00369 \
00370 if(ths->flags & STEEPEST_DESCENT) \
00371 fftw_free(ths->v_iter); \
00372 \
00373 fftw_free(ths->p_hat_iter); \
00374 fftw_free(ths->f_hat_iter); \
00375 \
00376 fftw_free(ths->r_iter); \
00377 fftw_free(ths->y); \
00378 }
00380 MACRO_SOLVER_IMPL(nfft, complex, double complex)
00381 MACRO_SOLVER_IMPL(nfct, double, double)
00382 MACRO_SOLVER_IMPL(nfst, double, double)
00383 MACRO_SOLVER_IMPL(nnfft, complex, double complex)
00384 MACRO_SOLVER_IMPL(mri_inh_2d1d, complex, double complex)
00385 MACRO_SOLVER_IMPL(mri_inh_3d, complex, double complex)
00386 MACRO_SOLVER_IMPL(nfsft, complex, double complex)
