Mercurial > octave-nkf

comparison liboctave/oct-fftw.cc @ 4775:88b638195bd1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2004-02-16 19:57:06 by jwe]
author jwe
date Mon, 16 Feb 2004 19:57:06 +0000
parents 0ff45249d321
children 55975a3073be
comparison
equal deleted inserted replaced
4774:0ff45249d321 4775:88b638195bd1
22 #include <config.h> 22 #include <config.h>
23 #endif 23 #endif
24 24
25 #if defined (HAVE_FFTW3) 25 #if defined (HAVE_FFTW3)
26 26
27 #include <iostream>
28 #include <vector>
29
27 #include "oct-fftw.h" 30 #include "oct-fftw.h"
28 #include "lo-error.h" 31 #include "lo-error.h"
29 #include <iostream>
30 32
31 // Helper class to create and cache fftw plans for both 1d and 2d. This 33 // Helper class to create and cache fftw plans for both 1d and 2d. This
32 // implementation uses FFTW_ESTIMATE to create the plans, which in theory 34 // implementation uses FFTW_ESTIMATE to create the plans, which in theory
33 // is suboptimal, but provides quite reasonable performance. 35 // is suboptimal, but provides quite reasonable performance.
34 36
105 const Complex *in, Complex *out) 107 const Complex *in, Complex *out)
106 { 108 {
107 int which = (dir == FFTW_FORWARD) ? 0 : 1; 109 int which = (dir == FFTW_FORWARD) ? 0 : 1;
108 fftw_plan *cur_plan_p = &plan[which]; 110 fftw_plan *cur_plan_p = &plan[which];
109 bool create_new_plan = false; 111 bool create_new_plan = false;
110 char in_align = (static_cast<int> (in)) & 0xF; 112 char in_align = (X_CAST (int, in)) & 0xF;
111 char out_align = (static_cast<int (out)) & 0xF; 113 char out_align = (X_CAST (int, out)) & 0xF;
112 114
113 if (plan[which] == 0 || d[which] != dist || s[which] != stride || 115 if (plan[which] == 0 || d[which] != dist || s[which] != stride ||
114 r[which] != rank || h[which] != howmany || 116 r[which] != rank || h[which] != howmany ||
115 ialign[which] != in_align || oalign[which] != out_align) 117 ialign[which] != in_align || oalign[which] != out_align)
116 create_new_plan = true; 118 create_new_plan = true;
159 int howmany, int stride, int dist, 161 int howmany, int stride, int dist,
160 const double *in, Complex *out) 162 const double *in, Complex *out)
161 { 163 {
162 fftw_plan *cur_plan_p = &rplan; 164 fftw_plan *cur_plan_p = &rplan;
163 bool create_new_plan = false; 165 bool create_new_plan = false;
164 char in_align = (static_cast<int> (in)) & 0xF; 166 char in_align = (X_CAST (int, in)) & 0xF;
165 char out_align = (static_cast<int (out)) & 0xF; 167 char out_align = (X_CAST (int, out)) & 0xF;
166 168
167 if (rplan == 0 || rd != dist || rs != stride || 169 if (rplan == 0 || rd != dist || rs != stride ||
168 rr != rank || rh != howmany || 170 rr != rank || rh != howmany ||
169 rialign != in_align || roalign != out_align) 171 rialign != in_align || roalign != out_align)
170 create_new_plan = true; 172 create_new_plan = true;
208 return *cur_plan_p; 210 return *cur_plan_p;
209 } 211 }
210 212
211 static octave_fftw_planner fftw_planner; 213 static octave_fftw_planner fftw_planner;
212 214
213 static inline void convert_packcomplex_1d (Complex *out, size_t nr, 215 static inline void
214 size_t nc, int stride, int dist) 216 convert_packcomplex_1d (Complex *out, size_t nr, size_t nc,
217 int stride, int dist)
215 { 218 {
216 // Fill in the missing data 219 // Fill in the missing data
217 for (size_t i = 0; i < nr; i++) 220 for (size_t i = 0; i < nr; i++)
218 for (size_t j = nc/2+1; j < nc; j++) 221 for (size_t j = nc/2+1; j < nc; j++)
219 out[j*stride + i*dist] = conj(out[(nc - j)*stride + i*dist]); 222 out[j*stride + i*dist] = conj(out[(nc - j)*stride + i*dist]);
220 } 223 }
221 224
222 static inline void convert_packcomplex_Nd (Complex *out, 225 static inline void
223 const dim_vector &dv) 226 convert_packcomplex_Nd (Complex *out, const dim_vector &dv)
224 { 227 {
225 size_t nc = dv(0); 228 size_t nc = dv(0);
226 size_t nr = dv(1); 229 size_t nr = dv(1);
227 size_t np = (dv.length() > 2 ? dv.numel () / nc / nr : 1); 230 size_t np = (dv.length() > 2 ? dv.numel () / nc / nr : 1);
228 size_t nrp = nr * np; 231 size_t nrp = nr * np;