NFFT  3.3.0
vector3.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: util.c 3483 2010-04-23 19:02:34Z keiner $ */
20 
21 #include "infft.h"
22 
24 void Y(upd_axpy_complex)(C *x, R a, C *y, INT n)
25 {
26  INT k;
27 
28  for (k = 0; k < n; k++)
29  x[k] = a * x[k] + y[k];
30 }
31 
33 void Y(upd_axpy_double)(R *x, R a, R *y, INT n)
34 {
35  INT k;
36 
37  for (k = 0; k < n; k++)
38  x[k] = a * x[k] + y[k];
39 }
40 
41 
43 void Y(upd_xpay_complex)(C *x, R a, C *y, INT n)
44 {
45  INT k;
46 
47  for (k = 0; k < n; k++)
48  x[k] += a * y[k];
49 }
50 
52 void Y(upd_xpay_double)(R *x, R a, R *y, INT n)
53 {
54  INT k;
55 
56  for (k = 0; k < n; k++)
57  x[k] += a * y[k];
58 }
59 
61 void Y(upd_axpby_complex)(C *x, R a, C *y, R b, INT n)
62 {
63  INT k;
64 
65  for (k = 0; k < n; k++)
66  x[k] = a * x[k] + b * y[k];
67 }
68 
70 void Y(upd_axpby_double)(R *x, R a, R *y, R b, INT n)
71 {
72  INT k;
73 
74  for (k = 0; k < n; k++)
75  x[k] = a * x[k] + b * y[k];
76 }
77 
79 void Y(upd_xpawy_complex)(C *x, R a, R *w, C *y, INT n)
80 {
81  INT k;
82 
83  for (k = 0; k < n; k++)
84  x[k] += a * w[k] * y[k];
85 }
86 
88 void Y(upd_xpawy_double)(R *x, R a, R *w, R *y, INT n)
89 {
90  INT k;
91 
92  for (k = 0; k < n; k++)
93  x[k] += a * w[k] * y[k];
94 }
95 
97 void Y(upd_axpwy_complex)(C *x, R a, R *w, C *y, INT n)
98 {
99  INT k;
100 
101  for (k = 0; k < n; k++)
102  x[k] = a * x[k] + w[k] * y[k];
103 }
104 
106 void Y(upd_axpwy_double)(R *x, R a, R *w, R *y, INT n)
107 {
108  INT k;
109 
110  for (k = 0; k < n; k++)
111  x[k] = a * x[k] + w[k] * y[k];
112 }
113 
115 void Y(fftshift_complex)(C *x, INT d, INT* N)
116 {
117  INT d_pre, d_act, d_post;
118  INT N_pre, N_act, N_post;
119  INT k_pre, k_act, k_post;
120  INT k, k_swap;
121 
122  C x_swap;
123 
124  for (d_act = 0; d_act < d; d_act++)
125  {
126  for (d_pre = 0, N_pre = 1; d_pre < d_act; d_pre++)
127  N_pre *= N[d_pre];
128 
129  N_act = N[d_act];
130 
131  for (d_post = d_act + 1, N_post = 1; d_post < d; d_post++)
132  N_post *= N[d_post];
133 
134  for (k_pre = 0; k_pre < N_pre; k_pre++)
135  for (k_act = 0; k_act < N_act / 2; k_act++)
136  for (k_post = 0; k_post < N_post; k_post++)
137  {
138  k = (k_pre * N_act + k_act) * N_post + k_post;
139  k_swap = (k_pre * N_act + k_act + N_act / 2) * N_post + k_post;
140 
141  x_swap = x[k];
142  x[k] = x[k_swap];
143  x[k_swap] = x_swap;
144  }
145  }
146 }
147 
149 void Y(fftshift_complex_int)(C *x, int d, int* N)
150 {
151  int d_pre, d_act, d_post;
152  int N_pre, N_act, N_post;
153  int k_pre, k_act, k_post;
154  int k, k_swap;
155 
156  C x_swap;
157 
158  for (d_act = 0; d_act < d; d_act++)
159  {
160  for (d_pre = 0, N_pre = 1; d_pre < d_act; d_pre++)
161  N_pre *= N[d_pre];
162 
163  N_act = N[d_act];
164 
165  for (d_post = d_act + 1, N_post = 1; d_post < d; d_post++)
166  N_post *= N[d_post];
167 
168  for (k_pre = 0; k_pre < N_pre; k_pre++)
169  for (k_act = 0; k_act < N_act / 2; k_act++)
170  for (k_post = 0; k_post < N_post; k_post++)
171  {
172  k = (k_pre * N_act + k_act) * N_post + k_post;
173  k_swap = (k_pre * N_act + k_act + N_act / 2) * N_post + k_post;
174 
175  x_swap = x[k];
176  x[k] = x[k_swap];
177  x[k_swap] = x_swap;
178  }
179  }
180 }