/src/math/nfft_f.c

https://gitlab.com/neelravi/octopus · C · 335 lines · 247 code · 64 blank · 24 comment · 17 complexity · 1067d461f44065a3cd29004689924d48 MD5 · raw file

  1. /*
  2. Copyright (C) 2002 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira
  3. This program is free software; you can redistribute it and/or modify
  4. it under the terms of the GNU General Public License as published by
  5. the Free Software Foundation; either version 2, or (at your option)
  6. any later version.
  7. This program is distributed in the hope that it will be useful,
  8. but WITHOUT ANY WARRANTY; without even the implied warranty of
  9. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  10. GNU General Public License for more details.
  11. You should have received a copy of the GNU General Public License
  12. along with this program; if not, write to the Free Software
  13. Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  14. 02110-1301, USA.
  15. $Id$
  16. */
  17. #include <config.h>
  18. #if defined(HAVE_NFFT)
  19. #include <stdio.h>
  20. #include <math.h>
  21. #include <stdlib.h>
  22. #include <string.h>
  23. #include <complex.h>
  24. #include "nfft3util.h"
  25. #include "nfft3.h"
  26. // NFFT FUNCTIONS
  27. void FC_FUNC(oct_nfft_init_1d,OCT_NFFT_INIT_1D)
  28. (nfft_plan *plan, int *N1, int *M)
  29. {
  30. nfft_init_1d(plan, *N1, *M);
  31. }
  32. void FC_FUNC(oct_nfft_init_2d,OCT_NFFT_INIT_2D)
  33. (nfft_plan *plan, int *N1, int *N2, int *M)
  34. {
  35. nfft_init_2d(plan, *N1, *N2, *M);
  36. }
  37. void FC_FUNC(oct_nfft_init_3d,OCT_NFFT_INIT_3D)
  38. (nfft_plan *plan, int *N1, int *N2, int *N3, int *M)
  39. {
  40. nfft_init_3d(plan, *N1, *N2, *N3, *M);
  41. }
  42. void FC_FUNC(oct_nfft_init_guru,OCT_NFFT_INIT_GURU)
  43. (nfft_plan *ths, int *d, int *N, int *M, int *n, int *m, unsigned *nfft_flags, unsigned *fftw_flags)
  44. {
  45. nfft_init_guru (ths, *d, N, *M, n, *m, *nfft_flags, *fftw_flags);
  46. }
  47. void FC_FUNC(oct_nfft_check, OCT_NFFT_CHECK)
  48. (nfft_plan *ths)
  49. {
  50. nfft_check (ths);
  51. }
  52. void FC_FUNC(oct_nfft_finalize,OCT_NFFT_FINALIZE)
  53. (nfft_plan *ths)
  54. {
  55. nfft_finalize (ths);
  56. }
  57. void FC_FUNC(oct_nfft_trafo,OCT_NFFT_TRAFO)
  58. (nfft_plan *ths)
  59. {
  60. nfft_trafo (ths);
  61. }
  62. void FC_FUNC(oct_nfft_adjoint,OCT_NFFT_ADJOINT)
  63. (nfft_plan *ths)
  64. {
  65. nfft_adjoint (ths);
  66. }
  67. void FC_FUNC(oct_nfft_precompute_one_psi_1d, OCT_NFFT_PRECOMPUTE_ONE_PSI_1D)
  68. (nfft_plan *plan, int *m, double *X1)
  69. {
  70. int ii;
  71. int M = *m;
  72. for (ii=0;ii< M;ii++){
  73. plan->x[ii] = X1[ii];
  74. }
  75. if(plan->nfft_flags & PRE_ONE_PSI)
  76. nfft_precompute_one_psi(plan);
  77. }
  78. void FC_FUNC(oct_nfft_precompute_one_psi_2d, OCT_NFFT_PRECOMPUTE_ONE_PSI_2D)
  79. (nfft_plan *plan, int *M, double* X1, double* X2)
  80. {
  81. int ii;
  82. int jj;
  83. for (ii=0; ii< M[0]; ii++){
  84. for (jj=0; jj< M[1]; jj++){
  85. plan->x[2*(M[1] * ii + jj) + 0] = X1[ii];
  86. plan->x[2*(M[1] * ii + jj) + 1] = X2[jj];
  87. }
  88. }
  89. if(plan->nfft_flags & PRE_ONE_PSI)
  90. nfft_precompute_one_psi(plan);
  91. }
  92. void FC_FUNC(oct_nfft_precompute_one_psi_3d, OCT_NFFT_PRECOMPUTE_ONE_PSI_3D)
  93. (nfft_plan *plan, int *M, double* X1, double* X2, double* X3)
  94. {
  95. int ii,jj,kk;
  96. for (ii=0;ii< M[0];ii++){
  97. for (jj=0;jj< M[1];jj++){
  98. for (kk=0;kk< M[2];kk++){
  99. plan->x[3*(M[1]*M[2]*ii + M[2]*jj + kk) + 0] = X1[ii];
  100. plan->x[3*(M[1]*M[2]*ii + M[2]*jj + kk) + 1] = X2[jj];
  101. plan->x[3*(M[1]*M[2]*ii + M[2]*jj + kk) + 2] = X3[kk];
  102. }
  103. }
  104. }
  105. if(plan->nfft_flags & PRE_ONE_PSI)
  106. nfft_precompute_one_psi(plan);
  107. }
  108. // Type dependent functions
  109. // ********** COMPLEX ************
  110. void FC_FUNC(zoct_set_f, ZOCT_SET_F)
  111. (nfft_plan *plan, int *M, int *DIM, double complex *VAL, int *IX, int *IY, int *IZ)
  112. {
  113. int dim = *DIM;
  114. int ix = *IX;
  115. int iy = *IY;
  116. int iz = *IZ;
  117. double complex val = *VAL;
  118. switch (dim){
  119. case 1:
  120. plan->f[ix-1] = val;
  121. break;
  122. case 2:
  123. plan->f[(ix-1)*M[1] + (iy-1)] = val;
  124. break;
  125. case 3:
  126. plan->f[(ix-1)*M[1]*M[2] + (iy-1)*M[2] + (iz-1)] = val;
  127. break;
  128. }
  129. }
  130. void FC_FUNC(zoct_get_f, ZOCT_GET_F)
  131. (nfft_plan *plan, int *M, int *DIM, double complex *val, int *IX, int *IY, int *IZ)
  132. {
  133. int dim = *DIM;
  134. int ix = *IX;
  135. int iy = *IY;
  136. int iz = *IZ;
  137. switch (dim){
  138. case 1:
  139. *val = plan->f[ix-1];
  140. break;
  141. case 2:
  142. *val = plan->f[(ix-1)*M[1] + (iy-1)];
  143. break;
  144. case 3:
  145. *val = plan->f[(ix-1)*M[1]*M[2] + (iy-1)*M[2] + (iz-1)];
  146. break;
  147. }
  148. }
  149. void FC_FUNC(zoct_set_f_hat, ZOCT_SET_F_HAT)
  150. (nfft_plan *plan, int *DIM, double complex *VAL, int *IX, int *IY, int *IZ)
  151. {
  152. int dim = *DIM;
  153. int ix = *IX;
  154. int iy = *IY;
  155. int iz = *IZ;
  156. double complex val = *VAL;
  157. switch (dim){
  158. case 1:
  159. plan->f_hat[ix-1] = val;
  160. break;
  161. case 2:
  162. plan->f_hat[(ix-1)*plan->N[1] + (iy-1)] = val;
  163. break;
  164. case 3:
  165. plan->f_hat[(ix-1)*plan->N[1]*plan->N[2] + (iy-1)*plan->N[2] + (iz-1)] = val;
  166. break;
  167. }
  168. }
  169. void FC_FUNC(zoct_get_f_hat, ZOCT_GET_F_HAT)
  170. (nfft_plan *plan, int *DIM, double complex *val, int *IX, int *IY, int *IZ)
  171. {
  172. int dim = *DIM;
  173. int ix = *IX;
  174. int iy = *IY;
  175. int iz = *IZ;
  176. switch (dim){
  177. case 1:
  178. *val = plan->f_hat[ix-1];
  179. break;
  180. case 2:
  181. *val = plan->f_hat[(ix-1)*plan->N[1] + (iy-1)];
  182. break;
  183. case 3:
  184. *val = plan->f_hat[(ix-1)*plan->N[1]*plan->N[2] + (iy-1)*plan->N[2] + (iz-1)];
  185. break;
  186. }
  187. }
  188. // ********** DOUBLE ************
  189. void FC_FUNC(doct_set_f, DOCT_SET_F)
  190. (nfft_plan *plan, int *M, int *DIM, double *VAL, int *IX, int *IY, int *IZ)
  191. {
  192. int dim = *DIM;
  193. int ix = *IX;
  194. int iy = *IY;
  195. int iz = *IZ;
  196. double val = *VAL;
  197. switch (dim){
  198. case 1:
  199. plan->f[ix-1] = val;
  200. break;
  201. case 2:
  202. plan->f[(ix-1)*M[1] + (iy-1)] = val;
  203. break;
  204. case 3:
  205. plan->f[(ix-1)*M[1]*M[2] + (iy-1)*M[2] + (iz-1)] = val;
  206. break;
  207. }
  208. }
  209. void FC_FUNC(doct_get_f, DOCT_GET_F)
  210. (nfft_plan *plan, int *M, int *DIM, double *val, int *IX, int *IY, int *IZ)
  211. {
  212. int dim = *DIM;
  213. int ix = *IX;
  214. int iy = *IY;
  215. int iz = *IZ;
  216. switch (dim){
  217. case 1:
  218. *val = plan->f[ix-1];
  219. break;
  220. case 2:
  221. *val = plan->f[(ix-1)*M[1] + (iy-1)];
  222. break;
  223. case 3:
  224. *val = plan->f[(ix-1)*M[1]*M[2] + (iy-1)*M[2] + (iz-1)];
  225. break;
  226. }
  227. }
  228. void FC_FUNC(doct_set_f_hat, DOCT_SET_F_HAT)
  229. (nfft_plan *plan, int *DIM, double *VAL, int *IX, int *IY, int *IZ)
  230. {
  231. int dim = *DIM;
  232. int ix = *IX;
  233. int iy = *IY;
  234. int iz = *IZ;
  235. double val = *VAL;
  236. switch (dim){
  237. case 1:
  238. plan->f_hat[ix-1] = val;
  239. break;
  240. case 2:
  241. plan->f_hat[(ix-1)*plan->N[1] + (iy-1)] = val;
  242. break;
  243. case 3:
  244. plan->f_hat[(ix-1)*plan->N[1]*plan->N[2] + (iy-1)*plan->N[2] + (iz-1)] = val;
  245. break;
  246. }
  247. }
  248. void FC_FUNC(doct_get_f_hat, DOCT_GET_F_HAT)
  249. (nfft_plan *plan, int *DIM, double *val, int *IX, int *IY, int *IZ)
  250. {
  251. int dim = *DIM;
  252. int ix = *IX;
  253. int iy = *IY;
  254. int iz = *IZ;
  255. switch (dim){
  256. case 1:
  257. *val = plan->f_hat[ix-1];
  258. break;
  259. case 2:
  260. *val = plan->f_hat[(ix-1)*plan->N[1] + (iy-1)];
  261. break;
  262. case 3:
  263. *val = plan->f_hat[(ix-1)*plan->N[1]*plan->N[2] + (iy-1)*plan->N[2] + (iz-1)];
  264. break;
  265. }
  266. }
  267. #endif