NFFT  3.3.1
vector3.c
00001 /*
00002  * Copyright (c) 2002, 2016 Jens Keiner, Stefan Kunis, Daniel Potts
00003  *
00004  * This program is free software; you can redistribute it and/or modify it under
00005  * the terms of the GNU General Public License as published by the Free Software
00006  * Foundation; either version 2 of the License, or (at your option) any later
00007  * version.
00008  *
00009  * This program is distributed in the hope that it will be useful, but WITHOUT
00010  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00011  * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
00012  * details.
00013  *
00014  * You should have received a copy of the GNU General Public License along with
00015  * this program; if not, write to the Free Software Foundation, Inc., 51
00016  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00017  */
00018 
00019 #include "infft.h"
00020 
00022 void Y(upd_axpy_complex)(C *x, R a, C *y, INT n)
00023 {
00024   INT k;
00025 
00026   for (k = 0; k < n; k++)
00027     x[k] = a * x[k] + y[k];
00028 }
00029 
00031 void Y(upd_axpy_double)(R *x, R a, R *y, INT n)
00032 {
00033   INT k;
00034 
00035   for (k = 0; k < n; k++)
00036     x[k] = a * x[k] + y[k];
00037 }
00038 
00039 
00041 void Y(upd_xpay_complex)(C *x, R a, C *y, INT n)
00042 {
00043   INT k;
00044 
00045   for (k = 0; k < n; k++)
00046     x[k] += a * y[k];
00047 }
00048 
00050 void Y(upd_xpay_double)(R *x, R a, R *y, INT n)
00051 {
00052   INT k;
00053 
00054   for (k = 0; k < n; k++)
00055     x[k] += a * y[k];
00056 }
00057 
00059 void Y(upd_axpby_complex)(C *x, R a, C *y, R b, INT n)
00060 {
00061   INT k;
00062 
00063   for (k = 0; k < n; k++)
00064     x[k] = a * x[k] + b * y[k];
00065 }
00066 
00068 void Y(upd_axpby_double)(R *x, R a, R *y, R b, INT n)
00069 {
00070   INT k;
00071 
00072   for (k = 0; k < n; k++)
00073     x[k] = a * x[k] + b * y[k];
00074 }
00075 
00077 void Y(upd_xpawy_complex)(C *x, R a, R *w, C *y, INT n)
00078 {
00079   INT k;
00080 
00081   for (k = 0; k < n; k++)
00082     x[k] += a * w[k] * y[k];
00083 }
00084 
00086 void Y(upd_xpawy_double)(R *x, R a, R *w, R *y, INT n)
00087 {
00088   INT k;
00089 
00090   for (k = 0; k < n; k++)
00091     x[k] += a * w[k] * y[k];
00092 }
00093 
00095 void Y(upd_axpwy_complex)(C *x, R a, R *w, C *y, INT n)
00096 {
00097   INT k;
00098 
00099   for (k = 0; k < n; k++)
00100     x[k] = a * x[k] + w[k] * y[k];
00101 }
00102 
00104 void Y(upd_axpwy_double)(R *x, R a, R *w, R *y, INT n)
00105 {
00106   INT k;
00107 
00108   for (k = 0; k < n; k++)
00109     x[k] = a * x[k] + w[k] * y[k];
00110 }
00111 
00113 void Y(fftshift_complex)(C *x, INT d, INT* N)
00114 {
00115   INT d_pre, d_act, d_post;
00116   INT N_pre, N_act, N_post;
00117   INT k_pre, k_act, k_post;
00118   INT k, k_swap;
00119 
00120   C x_swap;
00121 
00122   for (d_act = 0; d_act < d; d_act++)
00123   {
00124     for (d_pre = 0, N_pre = 1; d_pre < d_act; d_pre++)
00125       N_pre *= N[d_pre];
00126 
00127     N_act = N[d_act];
00128 
00129     for (d_post = d_act + 1, N_post = 1; d_post < d; d_post++)
00130       N_post *= N[d_post];
00131 
00132     for (k_pre = 0; k_pre < N_pre; k_pre++)
00133       for (k_act = 0; k_act < N_act / 2; k_act++)
00134         for (k_post = 0; k_post < N_post; k_post++)
00135         {
00136           k = (k_pre * N_act + k_act) * N_post + k_post;
00137           k_swap = (k_pre * N_act + k_act + N_act / 2) * N_post + k_post;
00138 
00139           x_swap = x[k];
00140           x[k] = x[k_swap];
00141           x[k_swap] = x_swap;
00142         }
00143   }
00144 }
00145 
00147 void Y(fftshift_complex_int)(C *x, int d, int* N)
00148 {
00149   int d_pre, d_act, d_post;
00150   int N_pre, N_act, N_post;
00151   int k_pre, k_act, k_post;
00152   int k, k_swap;
00153 
00154   C x_swap;
00155 
00156   for (d_act = 0; d_act < d; d_act++)
00157   {
00158     for (d_pre = 0, N_pre = 1; d_pre < d_act; d_pre++)
00159       N_pre *= N[d_pre];
00160 
00161     N_act = N[d_act];
00162 
00163     for (d_post = d_act + 1, N_post = 1; d_post < d; d_post++)
00164       N_post *= N[d_post];
00165 
00166     for (k_pre = 0; k_pre < N_pre; k_pre++)
00167       for (k_act = 0; k_act < N_act / 2; k_act++)
00168         for (k_post = 0; k_post < N_post; k_post++)
00169         {
00170           k = (k_pre * N_act + k_act) * N_post + k_post;
00171           k_swap = (k_pre * N_act + k_act + N_act / 2) * N_post + k_post;
00172 
00173           x_swap = x[k];
00174           x[k] = x[k_swap];
00175           x[k_swap] = x_swap;
00176         }
00177   }
00178 }