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