Mercurial > hg > octave-lyh

comparison liboctave/dNDArray.cc @ 4773:ccfbd6047a54

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2004-02-16 19:02:32 by jwe]
author jwe
date Mon, 16 Feb 2004 19:02:33 +0000
parents ef5e598f099b
children 5eb5b8aaed8a
comparison
equal deleted inserted replaced
4772:9eed17b2c8d1 4773:ccfbd6047a54
32 #include <cfloat> 32 #include <cfloat>
33 33
34 #include "Array-util.h" 34 #include "Array-util.h"
35 #include "dNDArray.h" 35 #include "dNDArray.h"
36 #include "mx-base.h" 36 #include "mx-base.h"
37 #include "f77-fcn.h"
37 #include "lo-error.h" 38 #include "lo-error.h"
38 #include "lo-ieee.h" 39 #include "lo-ieee.h"
39 #include "lo-mappers.h" 40 #include "lo-mappers.h"
41
42 #if defined (HAVE_FFTW3)
43 #include "oct-fftw.h"
44
45 ComplexNDArray
46 NDArray::fourier (int dim) const
47 {
48 dim_vector dv = dims ();
49
50 if (dim > dv.length () || dim < 0)
51 return ComplexNDArray ();
52
53 int stride = 1;
54 int n = dv(dim);
55
56 for (int i = 0; i < dim; i++)
57 stride *= dv(i);
58
59 int howmany = numel () / dv (dim);
60 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
61 int nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride);
62 int dist = (stride == 1 ? n : 1);
63
64 const double *in (fortran_vec ());
65 ComplexNDArray retval (dv);
66 Complex *out (retval.fortran_vec ());
67
68 // Need to be careful here about the distance between fft's
69 for (int k = 0; k < nloop; k++)
70 octave_fftw::fft (in + k * stride * n, out + k * stride * n,
71 n, howmany, stride, dist);
72
73 return retval;
74 }
75
76 ComplexNDArray
77 NDArray::ifourier (const int dim) const
78 {
79 dim_vector dv = dims ();
80
81 if (dim > dv.length () || dim < 0)
82 return ComplexNDArray ();
83
84 int stride = 1;
85 int n = dv(dim);
86
87 for (int i = 0; i < dim; i++)
88 stride *= dv(i);
89
90 int howmany = numel () / dv (dim);
91 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
92 int nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride);
93 int dist = (stride == 1 ? n : 1);
94
95 ComplexNDArray retval (*this);
96 Complex *out (retval.fortran_vec ());
97
98 // Need to be careful here about the distance between fft's
99 for (int k = 0; k < nloop; k++)
100 octave_fftw::ifft (out + k * stride * n, out + k * stride * n,
101 n, howmany, stride, dist);
102
103 return retval;
104 }
105
106 ComplexNDArray
107 NDArray::fourier2d (void) const
108 {
109 dim_vector dv = dims();
110 if (dv.length () < 2)
111 return ComplexNDArray ();
112
113 dim_vector dv2(dv(0), dv(1));
114 const double *in = fortran_vec ();
115 ComplexNDArray retval (dv);
116 Complex *out = retval.fortran_vec ();
117 int howmany = numel() / dv(0) / dv(1);
118 int dist = dv(0) * dv(1);
119
120 for (int i=0; i < howmany; i++)
121 octave_fftw::fftNd (in + i*dist, out + i*dist, 2, dv2);
122
123 return retval;
124 }
125
126 ComplexNDArray
127 NDArray::ifourier2d (void) const
128 {
129 dim_vector dv = dims();
130 if (dv.length () < 2)
131 return ComplexNDArray ();
132
133 dim_vector dv2(dv(0), dv(1));
134 ComplexNDArray retval (*this);
135 Complex *out = retval.fortran_vec ();
136 int howmany = numel() / dv(0) / dv(1);
137 int dist = dv(0) * dv(1);
138
139 for (int i=0; i < howmany; i++)
140 octave_fftw::ifftNd (out + i*dist, out + i*dist, 2, dv2);
141
142 return retval;
143 }
144
145 ComplexNDArray
146 NDArray::fourierNd (void) const
147 {
148 dim_vector dv = dims ();
149 int rank = dv.length ();
150
151 const double *in (fortran_vec ());
152 ComplexNDArray retval (dv);
153 Complex *out (retval.fortran_vec ());
154
155 octave_fftw::fftNd (in, out, rank, dv);
156
157 return retval;
158 }
159
160 ComplexNDArray
161 NDArray::ifourierNd (void) const
162 {
163 dim_vector dv = dims ();
164 int rank = dv.length ();
165
166 ComplexNDArray tmp (*this);
167 Complex *in (tmp.fortran_vec ());
168 ComplexNDArray retval (dv);
169 Complex *out (retval.fortran_vec ());
170
171 octave_fftw::ifftNd (in, out, rank, dv);
172
173 return retval;
174 }
175
176 #else
177
178 extern "C"
179 {
180 // Note that the original complex fft routines were not written for
181 // double complex arguments. They have been modified by adding an
182 // implicit double precision (a-h,o-z) statement at the beginning of
183 // each subroutine.
184
185 F77_RET_T
186 F77_FUNC (cffti, CFFTI) (const int&, Complex*);
187
188 F77_RET_T
189 F77_FUNC (cfftf, CFFTF) (const int&, Complex*, Complex*);
190
191 F77_RET_T
192 F77_FUNC (cfftb, CFFTB) (const int&, Complex*, Complex*);
193 }
194
195 ComplexNDArray
196 NDArray::fourier (int dim) const
197 {
198 dim_vector dv = dims ();
199
200 if (dim > dv.length () || dim < 0)
201 return ComplexNDArray ();
202
203 ComplexNDArray retval (dv);
204 int npts = dv(dim);
205 int nn = 4*npts+15;
206 Array<Complex> wsave (nn);
207 Complex *pwsave = wsave.fortran_vec ();
208
209 OCTAVE_LOCAL_BUFFER (Complex, tmp, npts);
210
211 int stride = 1;
212
213 for (int i = 0; i < dim; i++)
214 stride *= dv(i);
215
216 int howmany = numel () / npts;
217 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
218 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
219 int dist = (stride == 1 ? npts : 1);
220
221 F77_FUNC (cffti, CFFTI) (npts, pwsave);
222
223 for (int k = 0; k < nloop; k++)
224 {
225 for (int j = 0; j < howmany; j++)
226 {
227 OCTAVE_QUIT;
228
229 for (int i = 0; i < npts; i++)
230 tmp[i] = elem((i + k*npts)*stride + j*dist);
231
232 F77_FUNC (cfftf, CFFTF) (npts, tmp, pwsave);
233
234 for (int i = 0; i < npts; i++)
235 retval ((i + k*npts)*stride + j*dist) = tmp[i];
236 }
237 }
238
239 return retval;
240 }
241
242 ComplexNDArray
243 NDArray::ifourier (int dim) const
244 {
245 dim_vector dv = dims ();
246
247 if (dim > dv.length () || dim < 0)
248 return ComplexNDArray ();
249
250 ComplexNDArray retval (dv);
251 int npts = dv(dim);
252 int nn = 4*npts+15;
253 Array<Complex> wsave (nn);
254 Complex *pwsave = wsave.fortran_vec ();
255
256 OCTAVE_LOCAL_BUFFER (Complex, tmp, npts);
257
258 int stride = 1;
259
260 for (int i = 0; i < dim; i++)
261 stride *= dv(i);
262
263 int howmany = numel () / npts;
264 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
265 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
266 int dist = (stride == 1 ? npts : 1);
267
268 F77_FUNC (cffti, CFFTI) (npts, pwsave);
269
270 for (int k = 0; k < nloop; k++)
271 {
272 for (int j = 0; j < howmany; j++)
273 {
274 OCTAVE_QUIT;
275
276 for (int i = 0; i < npts; i++)
277 tmp[i] = elem((i + k*npts)*stride + j*dist);
278
279 F77_FUNC (cfftb, CFFTB) (npts, tmp, pwsave);
280
281 for (int i = 0; i < npts; i++)
282 retval ((i + k*npts)*stride + j*dist) = tmp[i] /
283 static_cast<double> (npts);
284 }
285 }
286
287 return retval;
288 }
289
290 ComplexNDArray
291 NDArray::fourier2d (void) const
292 {
293 dim_vector dv = dims();
294 dim_vector dv2 (dv(0), dv(1));
295 int rank = 2;
296 ComplexNDArray retval (*this);
297 int stride = 1;
298
299 for (int i = 0; i < rank; i++)
300 {
301 int npts = dv2(i);
302 int nn = 4*npts+15;
303 Array<Complex> wsave (nn);
304 Complex *pwsave = wsave.fortran_vec ();
305 Array<Complex> row (npts);
306 Complex *prow = row.fortran_vec ();
307
308 int howmany = numel () / npts;
309 howmany = (stride == 1 ? howmany :
310 (howmany > stride ? stride : howmany));
311 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
312 int dist = (stride == 1 ? npts : 1);
313
314 F77_FUNC (cffti, CFFTI) (npts, pwsave);
315
316 for (int k = 0; k < nloop; k++)
317 {
318 for (int j = 0; j < howmany; j++)
319 {
320 OCTAVE_QUIT;
321
322 for (int l = 0; l < npts; l++)
323 prow[l] = retval ((l + k*npts)*stride + j*dist);
324
325 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
326
327 for (int l = 0; l < npts; l++)
328 retval ((l + k*npts)*stride + j*dist) = prow[l];
329 }
330 }
331
332 stride *= dv2(i);
333 }
334
335 return retval;
336 }
337
338 ComplexNDArray
339 NDArray::ifourier2d (void) const
340 {
341 dim_vector dv = dims();
342 dim_vector dv2 (dv(0), dv(1));
343 int rank = 2;
344 ComplexNDArray retval (*this);
345 int stride = 1;
346
347 for (int i = 0; i < rank; i++)
348 {
349 int npts = dv2(i);
350 int nn = 4*npts+15;
351 Array<Complex> wsave (nn);
352 Complex *pwsave = wsave.fortran_vec ();
353 Array<Complex> row (npts);
354 Complex *prow = row.fortran_vec ();
355
356 int howmany = numel () / npts;
357 howmany = (stride == 1 ? howmany :
358 (howmany > stride ? stride : howmany));
359 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
360 int dist = (stride == 1 ? npts : 1);
361
362 F77_FUNC (cffti, CFFTI) (npts, pwsave);
363
364 for (int k = 0; k < nloop; k++)
365 {
366 for (int j = 0; j < howmany; j++)
367 {
368 OCTAVE_QUIT;
369
370 for (int l = 0; l < npts; l++)
371 prow[l] = retval ((l + k*npts)*stride + j*dist);
372
373 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
374
375 for (int l = 0; l < npts; l++)
376 retval ((l + k*npts)*stride + j*dist) = prow[l] /
377 static_cast<double> (npts);
378 }
379 }
380
381 stride *= dv2(i);
382 }
383
384 return retval;
385 }
386
387 ComplexNDArray
388 NDArray::fourierNd (void) const
389 {
390 dim_vector dv = dims ();
391 int rank = dv.length ();
392 ComplexNDArray retval (*this);
393 int stride = 1;
394
395 for (int i = 0; i < rank; i++)
396 {
397 int npts = dv(i);
398 int nn = 4*npts+15;
399 Array<Complex> wsave (nn);
400 Complex *pwsave = wsave.fortran_vec ();
401 Array<Complex> row (npts);
402 Complex *prow = row.fortran_vec ();
403
404 int howmany = numel () / npts;
405 howmany = (stride == 1 ? howmany :
406 (howmany > stride ? stride : howmany));
407 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
408 int dist = (stride == 1 ? npts : 1);
409
410 F77_FUNC (cffti, CFFTI) (npts, pwsave);
411
412 for (int k = 0; k < nloop; k++)
413 {
414 for (int j = 0; j < howmany; j++)
415 {
416 OCTAVE_QUIT;
417
418 for (int l = 0; l < npts; l++)
419 prow[l] = retval ((l + k*npts)*stride + j*dist);
420
421 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
422
423 for (int l = 0; l < npts; l++)
424 retval ((l + k*npts)*stride + j*dist) = prow[l];
425 }
426 }
427
428 stride *= dv(i);
429 }
430
431 return retval;
432 }
433
434 ComplexNDArray
435 NDArray::ifourierNd (void) const
436 {
437 dim_vector dv = dims ();
438 int rank = dv.length ();
439 ComplexNDArray retval (*this);
440 int stride = 1;
441
442 for (int i = 0; i < rank; i++)
443 {
444 int npts = dv(i);
445 int nn = 4*npts+15;
446 Array<Complex> wsave (nn);
447 Complex *pwsave = wsave.fortran_vec ();
448 Array<Complex> row (npts);
449 Complex *prow = row.fortran_vec ();
450
451 int howmany = numel () / npts;
452 howmany = (stride == 1 ? howmany :
453 (howmany > stride ? stride : howmany));
454 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
455 int dist = (stride == 1 ? npts : 1);
456
457 F77_FUNC (cffti, CFFTI) (npts, pwsave);
458
459 for (int k = 0; k < nloop; k++)
460 {
461 for (int j = 0; j < howmany; j++)
462 {
463 OCTAVE_QUIT;
464
465 for (int l = 0; l < npts; l++)
466 prow[l] = retval ((l + k*npts)*stride + j*dist);
467
468 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
469
470 for (int l = 0; l < npts; l++)
471 retval ((l + k*npts)*stride + j*dist) = prow[l] /
472 static_cast<double> (npts);
473 }
474 }
475
476 stride *= dv(i);
477 }
478
479 return retval;
480 }
481
482 #endif
40 483
41 NDArray::NDArray (const boolNDArray& a) 484 NDArray::NDArray (const boolNDArray& a)
42 : MArrayN<double> (a.dims ()) 485 : MArrayN<double> (a.dims ())
43 { 486 {
44 for (int i = 0; i < a.length (); i++) 487 for (int i = 0; i < a.length (); i++)