|
4514
|
1 // N-D Array manipulations. |
|
|
2 /* |
|
|
3 |
|
|
4 Copyright (C) 1996, 1997 John W. Eaton |
|
|
5 |
|
|
6 This file is part of Octave. |
|
|
7 |
|
|
8 Octave is free software; you can redistribute it and/or modify it |
|
|
9 under the terms of the GNU General Public License as published by the |
|
|
10 Free Software Foundation; either version 2, or (at your option) any |
|
|
11 later version. |
|
|
12 |
|
|
13 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
16 for more details. |
|
|
17 |
|
|
18 You should have received a copy of the GNU General Public License |
|
|
19 along with Octave; see the file COPYING. If not, write to the Free |
|
|
20 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
|
|
21 |
|
|
22 */ |
|
|
23 |
|
|
24 #ifdef HAVE_CONFIG_H |
|
|
25 #include <config.h> |
|
|
26 #endif |
|
|
27 |
|
4687
|
28 #include <cfloat> |
|
4780
|
29 #include <vector> |
|
4687
|
30 |
|
4588
|
31 #include "Array-util.h" |
|
4514
|
32 #include "CNDArray.h" |
|
|
33 #include "mx-base.h" |
|
4773
|
34 #include "f77-fcn.h" |
|
4514
|
35 #include "lo-ieee.h" |
|
4687
|
36 #include "lo-mappers.h" |
|
4514
|
37 |
|
4773
|
38 #if defined (HAVE_FFTW3) |
|
4780
|
39 #include "oct-fftw.h" |
|
4773
|
40 #else |
|
|
41 extern "C" |
|
|
42 { |
|
|
43 // Note that the original complex fft routines were not written for |
|
|
44 // double complex arguments. They have been modified by adding an |
|
|
45 // implicit double precision (a-h,o-z) statement at the beginning of |
|
|
46 // each subroutine. |
|
|
47 |
|
|
48 F77_RET_T |
|
|
49 F77_FUNC (cffti, CFFTI) (const int&, Complex*); |
|
|
50 |
|
|
51 F77_RET_T |
|
|
52 F77_FUNC (cfftf, CFFTF) (const int&, Complex*, Complex*); |
|
|
53 |
|
|
54 F77_RET_T |
|
|
55 F77_FUNC (cfftb, CFFTB) (const int&, Complex*, Complex*); |
|
|
56 } |
|
|
57 #endif |
|
|
58 |
|
|
59 #if defined (HAVE_FFTW3) |
|
|
60 ComplexNDArray |
|
|
61 ComplexNDArray::fourier (int dim) const |
|
|
62 { |
|
|
63 dim_vector dv = dims (); |
|
|
64 |
|
|
65 if (dim > dv.length () || dim < 0) |
|
|
66 return ComplexNDArray (); |
|
|
67 |
|
|
68 int stride = 1; |
|
|
69 int n = dv(dim); |
|
|
70 |
|
|
71 for (int i = 0; i < dim; i++) |
|
|
72 stride *= dv(i); |
|
|
73 |
|
|
74 int howmany = numel () / dv (dim); |
|
|
75 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany)); |
|
|
76 int nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride); |
|
|
77 int dist = (stride == 1 ? n : 1); |
|
|
78 |
|
|
79 const Complex *in (fortran_vec ()); |
|
|
80 ComplexNDArray retval (dv); |
|
|
81 Complex *out (retval.fortran_vec ()); |
|
|
82 |
|
|
83 // Need to be careful here about the distance between fft's |
|
|
84 for (int k = 0; k < nloop; k++) |
|
|
85 octave_fftw::fft (in + k * stride * n, out + k * stride * n, |
|
|
86 n, howmany, stride, dist); |
|
|
87 |
|
|
88 return retval; |
|
|
89 } |
|
|
90 |
|
|
91 ComplexNDArray |
|
4816
|
92 ComplexNDArray::ifourier (int dim) const |
|
4773
|
93 { |
|
|
94 dim_vector dv = dims (); |
|
|
95 |
|
|
96 if (dim > dv.length () || dim < 0) |
|
|
97 return ComplexNDArray (); |
|
|
98 |
|
|
99 int stride = 1; |
|
|
100 int n = dv(dim); |
|
|
101 |
|
|
102 for (int i = 0; i < dim; i++) |
|
|
103 stride *= dv(i); |
|
|
104 |
|
|
105 int howmany = numel () / dv (dim); |
|
|
106 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany)); |
|
|
107 int nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride); |
|
|
108 int dist = (stride == 1 ? n : 1); |
|
|
109 |
|
|
110 const Complex *in (fortran_vec ()); |
|
|
111 ComplexNDArray retval (dv); |
|
|
112 Complex *out (retval.fortran_vec ()); |
|
|
113 |
|
|
114 // Need to be careful here about the distance between fft's |
|
|
115 for (int k = 0; k < nloop; k++) |
|
|
116 octave_fftw::ifft (in + k * stride * n, out + k * stride * n, |
|
|
117 n, howmany, stride, dist); |
|
|
118 |
|
|
119 return retval; |
|
|
120 } |
|
|
121 |
|
|
122 ComplexNDArray |
|
|
123 ComplexNDArray::fourier2d (void) const |
|
|
124 { |
|
|
125 dim_vector dv = dims(); |
|
|
126 if (dv.length () < 2) |
|
|
127 return ComplexNDArray (); |
|
|
128 |
|
|
129 dim_vector dv2(dv(0), dv(1)); |
|
|
130 const Complex *in = fortran_vec (); |
|
|
131 ComplexNDArray retval (dv); |
|
|
132 Complex *out = retval.fortran_vec (); |
|
|
133 int howmany = numel() / dv(0) / dv(1); |
|
|
134 int dist = dv(0) * dv(1); |
|
|
135 |
|
|
136 for (int i=0; i < howmany; i++) |
|
|
137 octave_fftw::fftNd (in + i*dist, out + i*dist, 2, dv2); |
|
|
138 |
|
|
139 return retval; |
|
|
140 } |
|
|
141 |
|
|
142 ComplexNDArray |
|
|
143 ComplexNDArray::ifourier2d (void) const |
|
|
144 { |
|
|
145 dim_vector dv = dims(); |
|
|
146 if (dv.length () < 2) |
|
|
147 return ComplexNDArray (); |
|
|
148 |
|
|
149 dim_vector dv2(dv(0), dv(1)); |
|
|
150 const Complex *in = fortran_vec (); |
|
|
151 ComplexNDArray retval (dv); |
|
|
152 Complex *out = retval.fortran_vec (); |
|
|
153 int howmany = numel() / dv(0) / dv(1); |
|
|
154 int dist = dv(0) * dv(1); |
|
|
155 |
|
|
156 for (int i=0; i < howmany; i++) |
|
|
157 octave_fftw::ifftNd (in + i*dist, out + i*dist, 2, dv2); |
|
|
158 |
|
|
159 return retval; |
|
|
160 } |
|
|
161 |
|
|
162 ComplexNDArray |
|
|
163 ComplexNDArray::fourierNd (void) const |
|
|
164 { |
|
|
165 dim_vector dv = dims (); |
|
|
166 int rank = dv.length (); |
|
|
167 |
|
|
168 const Complex *in (fortran_vec ()); |
|
|
169 ComplexNDArray retval (dv); |
|
|
170 Complex *out (retval.fortran_vec ()); |
|
|
171 |
|
|
172 octave_fftw::fftNd (in, out, rank, dv); |
|
|
173 |
|
|
174 return retval; |
|
|
175 } |
|
|
176 |
|
|
177 ComplexNDArray |
|
|
178 ComplexNDArray::ifourierNd (void) const |
|
|
179 { |
|
|
180 dim_vector dv = dims (); |
|
|
181 int rank = dv.length (); |
|
|
182 |
|
|
183 const Complex *in (fortran_vec ()); |
|
|
184 ComplexNDArray retval (dv); |
|
|
185 Complex *out (retval.fortran_vec ()); |
|
|
186 |
|
|
187 octave_fftw::ifftNd (in, out, rank, dv); |
|
|
188 |
|
|
189 return retval; |
|
|
190 } |
|
|
191 |
|
|
192 #else |
|
|
193 ComplexNDArray |
|
|
194 ComplexNDArray::fourier (int dim) const |
|
|
195 { |
|
|
196 dim_vector dv = dims (); |
|
|
197 |
|
|
198 if (dim > dv.length () || dim < 0) |
|
|
199 return ComplexNDArray (); |
|
|
200 |
|
|
201 ComplexNDArray retval (dv); |
|
|
202 int npts = dv(dim); |
|
|
203 int nn = 4*npts+15; |
|
|
204 Array<Complex> wsave (nn); |
|
|
205 Complex *pwsave = wsave.fortran_vec (); |
|
|
206 |
|
|
207 OCTAVE_LOCAL_BUFFER (Complex, tmp, npts); |
|
|
208 |
|
|
209 int stride = 1; |
|
|
210 |
|
|
211 for (int i = 0; i < dim; i++) |
|
|
212 stride *= dv(i); |
|
|
213 |
|
|
214 int howmany = numel () / npts; |
|
|
215 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany)); |
|
|
216 int nloop = (stride == 1 ? 1 : numel () / npts / stride); |
|
|
217 int dist = (stride == 1 ? npts : 1); |
|
|
218 |
|
|
219 F77_FUNC (cffti, CFFTI) (npts, pwsave); |
|
|
220 |
|
|
221 for (int k = 0; k < nloop; k++) |
|
|
222 { |
|
|
223 for (int j = 0; j < howmany; j++) |
|
|
224 { |
|
|
225 OCTAVE_QUIT; |
|
|
226 |
|
|
227 for (int i = 0; i < npts; i++) |
|
|
228 tmp[i] = elem((i + k*npts)*stride + j*dist); |
|
|
229 |
|
|
230 F77_FUNC (cfftf, CFFTF) (npts, tmp, pwsave); |
|
|
231 |
|
|
232 for (int i = 0; i < npts; i++) |
|
|
233 retval ((i + k*npts)*stride + j*dist) = tmp[i]; |
|
|
234 } |
|
|
235 } |
|
|
236 |
|
|
237 return retval; |
|
|
238 } |
|
|
239 |
|
|
240 ComplexNDArray |
|
|
241 ComplexNDArray::ifourier (int dim) const |
|
|
242 { |
|
|
243 dim_vector dv = dims (); |
|
|
244 |
|
|
245 if (dim > dv.length () || dim < 0) |
|
|
246 return ComplexNDArray (); |
|
|
247 |
|
|
248 ComplexNDArray retval (dv); |
|
|
249 int npts = dv(dim); |
|
|
250 int nn = 4*npts+15; |
|
|
251 Array<Complex> wsave (nn); |
|
|
252 Complex *pwsave = wsave.fortran_vec (); |
|
|
253 |
|
|
254 OCTAVE_LOCAL_BUFFER (Complex, tmp, npts); |
|
|
255 |
|
|
256 int stride = 1; |
|
|
257 |
|
|
258 for (int i = 0; i < dim; i++) |
|
|
259 stride *= dv(i); |
|
|
260 |
|
|
261 int howmany = numel () / npts; |
|
|
262 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany)); |
|
|
263 int nloop = (stride == 1 ? 1 : numel () / npts / stride); |
|
|
264 int dist = (stride == 1 ? npts : 1); |
|
|
265 |
|
|
266 F77_FUNC (cffti, CFFTI) (npts, pwsave); |
|
|
267 |
|
|
268 for (int k = 0; k < nloop; k++) |
|
|
269 { |
|
|
270 for (int j = 0; j < howmany; j++) |
|
|
271 { |
|
|
272 OCTAVE_QUIT; |
|
|
273 |
|
|
274 for (int i = 0; i < npts; i++) |
|
|
275 tmp[i] = elem((i + k*npts)*stride + j*dist); |
|
|
276 |
|
|
277 F77_FUNC (cfftb, CFFTB) (npts, tmp, pwsave); |
|
|
278 |
|
|
279 for (int i = 0; i < npts; i++) |
|
|
280 retval ((i + k*npts)*stride + j*dist) = tmp[i] / |
|
|
281 static_cast<double> (npts); |
|
|
282 } |
|
|
283 } |
|
|
284 |
|
|
285 return retval; |
|
|
286 } |
|
|
287 |
|
|
288 ComplexNDArray |
|
|
289 ComplexNDArray::fourier2d (void) const |
|
|
290 { |
|
|
291 dim_vector dv = dims (); |
|
|
292 dim_vector dv2 (dv(0), dv(1)); |
|
|
293 int rank = 2; |
|
|
294 ComplexNDArray retval (*this); |
|
|
295 int stride = 1; |
|
|
296 |
|
|
297 for (int i = 0; i < rank; i++) |
|
|
298 { |
|
|
299 int npts = dv2(i); |
|
|
300 int nn = 4*npts+15; |
|
|
301 Array<Complex> wsave (nn); |
|
|
302 Complex *pwsave = wsave.fortran_vec (); |
|
|
303 Array<Complex> row (npts); |
|
|
304 Complex *prow = row.fortran_vec (); |
|
|
305 |
|
|
306 int howmany = numel () / npts; |
|
|
307 howmany = (stride == 1 ? howmany : |
|
|
308 (howmany > stride ? stride : howmany)); |
|
|
309 int nloop = (stride == 1 ? 1 : numel () / npts / stride); |
|
|
310 int dist = (stride == 1 ? npts : 1); |
|
|
311 |
|
|
312 F77_FUNC (cffti, CFFTI) (npts, pwsave); |
|
|
313 |
|
|
314 for (int k = 0; k < nloop; k++) |
|
|
315 { |
|
|
316 for (int j = 0; j < howmany; j++) |
|
|
317 { |
|
|
318 OCTAVE_QUIT; |
|
|
319 |
|
|
320 for (int l = 0; l < npts; l++) |
|
|
321 prow[l] = retval ((l + k*npts)*stride + j*dist); |
|
|
322 |
|
|
323 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave); |
|
|
324 |
|
|
325 for (int l = 0; l < npts; l++) |
|
|
326 retval ((l + k*npts)*stride + j*dist) = prow[l]; |
|
|
327 } |
|
|
328 } |
|
|
329 |
|
|
330 stride *= dv2(i); |
|
|
331 } |
|
|
332 |
|
|
333 return retval; |
|
|
334 } |
|
|
335 |
|
|
336 ComplexNDArray |
|
|
337 ComplexNDArray::ifourier2d (void) const |
|
|
338 { |
|
|
339 dim_vector dv = dims(); |
|
|
340 dim_vector dv2 (dv(0), dv(1)); |
|
|
341 int rank = 2; |
|
|
342 ComplexNDArray retval (*this); |
|
|
343 int stride = 1; |
|
|
344 |
|
|
345 for (int i = 0; i < rank; i++) |
|
|
346 { |
|
|
347 int npts = dv2(i); |
|
|
348 int nn = 4*npts+15; |
|
|
349 Array<Complex> wsave (nn); |
|
|
350 Complex *pwsave = wsave.fortran_vec (); |
|
|
351 Array<Complex> row (npts); |
|
|
352 Complex *prow = row.fortran_vec (); |
|
|
353 |
|
|
354 int howmany = numel () / npts; |
|
|
355 howmany = (stride == 1 ? howmany : |
|
|
356 (howmany > stride ? stride : howmany)); |
|
|
357 int nloop = (stride == 1 ? 1 : numel () / npts / stride); |
|
|
358 int dist = (stride == 1 ? npts : 1); |
|
|
359 |
|
|
360 F77_FUNC (cffti, CFFTI) (npts, pwsave); |
|
|
361 |
|
|
362 for (int k = 0; k < nloop; k++) |
|
|
363 { |
|
|
364 for (int j = 0; j < howmany; j++) |
|
|
365 { |
|
|
366 OCTAVE_QUIT; |
|
|
367 |
|
|
368 for (int l = 0; l < npts; l++) |
|
|
369 prow[l] = retval ((l + k*npts)*stride + j*dist); |
|
|
370 |
|
|
371 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave); |
|
|
372 |
|
|
373 for (int l = 0; l < npts; l++) |
|
|
374 retval ((l + k*npts)*stride + j*dist) = prow[l] / |
|
|
375 static_cast<double> (npts); |
|
|
376 } |
|
|
377 } |
|
|
378 |
|
|
379 stride *= dv2(i); |
|
|
380 } |
|
|
381 |
|
|
382 return retval; |
|
|
383 } |
|
|
384 |
|
|
385 ComplexNDArray |
|
|
386 ComplexNDArray::fourierNd (void) const |
|
|
387 { |
|
|
388 dim_vector dv = dims (); |
|
|
389 int rank = dv.length (); |
|
|
390 ComplexNDArray retval (*this); |
|
|
391 int stride = 1; |
|
|
392 |
|
|
393 for (int i = 0; i < rank; i++) |
|
|
394 { |
|
|
395 int npts = dv(i); |
|
|
396 int nn = 4*npts+15; |
|
|
397 Array<Complex> wsave (nn); |
|
|
398 Complex *pwsave = wsave.fortran_vec (); |
|
|
399 Array<Complex> row (npts); |
|
|
400 Complex *prow = row.fortran_vec (); |
|
|
401 |
|
|
402 int howmany = numel () / npts; |
|
|
403 howmany = (stride == 1 ? howmany : |
|
|
404 (howmany > stride ? stride : howmany)); |
|
|
405 int nloop = (stride == 1 ? 1 : numel () / npts / stride); |
|
|
406 int dist = (stride == 1 ? npts : 1); |
|
|
407 |
|
|
408 F77_FUNC (cffti, CFFTI) (npts, pwsave); |
|
|
409 |
|
|
410 for (int k = 0; k < nloop; k++) |
|
|
411 { |
|
|
412 for (int j = 0; j < howmany; j++) |
|
|
413 { |
|
|
414 OCTAVE_QUIT; |
|
|
415 |
|
|
416 for (int l = 0; l < npts; l++) |
|
|
417 prow[l] = retval ((l + k*npts)*stride + j*dist); |
|
|
418 |
|
|
419 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave); |
|
|
420 |
|
|
421 for (int l = 0; l < npts; l++) |
|
|
422 retval ((l + k*npts)*stride + j*dist) = prow[l]; |
|
|
423 } |
|
|
424 } |
|
|
425 |
|
|
426 stride *= dv(i); |
|
|
427 } |
|
|
428 |
|
|
429 return retval; |
|
|
430 } |
|
|
431 |
|
|
432 ComplexNDArray |
|
|
433 ComplexNDArray::ifourierNd (void) const |
|
|
434 { |
|
|
435 dim_vector dv = dims (); |
|
|
436 int rank = dv.length (); |
|
|
437 ComplexNDArray retval (*this); |
|
|
438 int stride = 1; |
|
|
439 |
|
|
440 for (int i = 0; i < rank; i++) |
|
|
441 { |
|
|
442 int npts = dv(i); |
|
|
443 int nn = 4*npts+15; |
|
|
444 Array<Complex> wsave (nn); |
|
|
445 Complex *pwsave = wsave.fortran_vec (); |
|
|
446 Array<Complex> row (npts); |
|
|
447 Complex *prow = row.fortran_vec (); |
|
|
448 |
|
|
449 int howmany = numel () / npts; |
|
|
450 howmany = (stride == 1 ? howmany : |
|
|
451 (howmany > stride ? stride : howmany)); |
|
|
452 int nloop = (stride == 1 ? 1 : numel () / npts / stride); |
|
|
453 int dist = (stride == 1 ? npts : 1); |
|
|
454 |
|
|
455 F77_FUNC (cffti, CFFTI) (npts, pwsave); |
|
|
456 |
|
|
457 for (int k = 0; k < nloop; k++) |
|
|
458 { |
|
|
459 for (int j = 0; j < howmany; j++) |
|
|
460 { |
|
|
461 OCTAVE_QUIT; |
|
|
462 |
|
|
463 for (int l = 0; l < npts; l++) |
|
|
464 prow[l] = retval ((l + k*npts)*stride + j*dist); |
|
|
465 |
|
|
466 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave); |
|
|
467 |
|
|
468 for (int l = 0; l < npts; l++) |
|
|
469 retval ((l + k*npts)*stride + j*dist) = prow[l] / |
|
|
470 static_cast<double> (npts); |
|
|
471 } |
|
|
472 } |
|
|
473 |
|
|
474 stride *= dv(i); |
|
|
475 } |
|
|
476 |
|
|
477 return retval; |
|
|
478 } |
|
|
479 |
|
|
480 #endif |
|
|
481 |
|
4543
|
482 // unary operations |
|
|
483 |
|
|
484 boolNDArray |
|
|
485 ComplexNDArray::operator ! (void) const |
|
|
486 { |
|
|
487 boolNDArray b (dims ()); |
|
|
488 |
|
|
489 for (int i = 0; i < length (); i++) |
|
5139
|
490 b.elem (i) = elem (i) == 0.0; |
|
4543
|
491 |
|
|
492 return b; |
|
|
493 } |
|
|
494 |
|
4514
|
495 // XXX FIXME XXX -- this is not quite the right thing. |
|
|
496 |
|
4687
|
497 bool |
|
|
498 ComplexNDArray::any_element_is_inf_or_nan (void) const |
|
|
499 { |
|
|
500 int nel = nelem (); |
|
|
501 |
|
|
502 for (int i = 0; i < nel; i++) |
|
|
503 { |
|
|
504 Complex val = elem (i); |
|
|
505 if (xisinf (val) || xisnan (val)) |
|
|
506 return true; |
|
|
507 } |
|
|
508 return false; |
|
|
509 } |
|
|
510 |
|
|
511 // Return true if no elements have imaginary components. |
|
|
512 |
|
|
513 bool |
|
|
514 ComplexNDArray::all_elements_are_real (void) const |
|
|
515 { |
|
|
516 int nel = nelem (); |
|
|
517 |
|
|
518 for (int i = 0; i < nel; i++) |
|
|
519 { |
|
|
520 double ip = imag (elem (i)); |
|
|
521 |
|
|
522 if (ip != 0.0 || lo_ieee_signbit (ip)) |
|
|
523 return false; |
|
|
524 } |
|
|
525 |
|
|
526 return true; |
|
|
527 } |
|
|
528 |
|
|
529 // Return nonzero if any element of CM has a non-integer real or |
|
|
530 // imaginary part. Also extract the largest and smallest (real or |
|
|
531 // imaginary) values and return them in MAX_VAL and MIN_VAL. |
|
|
532 |
|
|
533 bool |
|
|
534 ComplexNDArray::all_integers (double& max_val, double& min_val) const |
|
|
535 { |
|
|
536 int nel = nelem (); |
|
|
537 |
|
|
538 if (nel > 0) |
|
|
539 { |
|
|
540 Complex val = elem (0); |
|
|
541 |
|
|
542 double r_val = real (val); |
|
|
543 double i_val = imag (val); |
|
|
544 |
|
|
545 max_val = r_val; |
|
|
546 min_val = r_val; |
|
|
547 |
|
|
548 if (i_val > max_val) |
|
|
549 max_val = i_val; |
|
|
550 |
|
|
551 if (i_val < max_val) |
|
|
552 min_val = i_val; |
|
|
553 } |
|
|
554 else |
|
|
555 return false; |
|
|
556 |
|
|
557 for (int i = 0; i < nel; i++) |
|
|
558 { |
|
|
559 Complex val = elem (i); |
|
|
560 |
|
|
561 double r_val = real (val); |
|
|
562 double i_val = imag (val); |
|
|
563 |
|
|
564 if (r_val > max_val) |
|
|
565 max_val = r_val; |
|
|
566 |
|
|
567 if (i_val > max_val) |
|
|
568 max_val = i_val; |
|
|
569 |
|
|
570 if (r_val < min_val) |
|
|
571 min_val = r_val; |
|
|
572 |
|
|
573 if (i_val < min_val) |
|
|
574 min_val = i_val; |
|
|
575 |
|
|
576 if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val) |
|
|
577 return false; |
|
|
578 } |
|
|
579 |
|
|
580 return true; |
|
|
581 } |
|
|
582 |
|
|
583 bool |
|
|
584 ComplexNDArray::too_large_for_float (void) const |
|
|
585 { |
|
|
586 int nel = nelem (); |
|
|
587 |
|
|
588 for (int i = 0; i < nel; i++) |
|
|
589 { |
|
|
590 Complex val = elem (i); |
|
|
591 |
|
|
592 double r_val = real (val); |
|
|
593 double i_val = imag (val); |
|
|
594 |
|
|
595 if (r_val > FLT_MAX |
|
|
596 || i_val > FLT_MAX |
|
|
597 || r_val < FLT_MIN |
|
|
598 || i_val < FLT_MIN) |
|
|
599 return true; |
|
|
600 } |
|
|
601 |
|
|
602 return false; |
|
|
603 } |
|
|
604 |
|
4556
|
605 boolNDArray |
|
4514
|
606 ComplexNDArray::all (int dim) const |
|
|
607 { |
|
4569
|
608 MX_ND_ANY_ALL_REDUCTION |
|
|
609 (MX_ND_ALL_EVAL (elem (iter_idx) == Complex (0, 0)), true); |
|
4514
|
610 } |
|
|
611 |
|
4556
|
612 boolNDArray |
|
4514
|
613 ComplexNDArray::any (int dim) const |
|
|
614 { |
|
4569
|
615 MX_ND_ANY_ALL_REDUCTION |
|
5110
|
616 (MX_ND_ANY_EVAL (elem (iter_idx) != Complex (0, 0) |
|
|
617 && ! (lo_ieee_isnan (::real (elem (iter_idx))) |
|
|
618 || lo_ieee_isnan (::imag (elem (iter_idx))))), |
|
|
619 false); |
|
4569
|
620 } |
|
|
621 |
|
4584
|
622 ComplexNDArray |
|
4569
|
623 ComplexNDArray::cumprod (int dim) const |
|
|
624 { |
|
4584
|
625 MX_ND_CUMULATIVE_OP (ComplexNDArray, Complex, Complex (1, 0), *); |
|
4569
|
626 } |
|
|
627 |
|
4584
|
628 ComplexNDArray |
|
4569
|
629 ComplexNDArray::cumsum (int dim) const |
|
|
630 { |
|
4584
|
631 MX_ND_CUMULATIVE_OP (ComplexNDArray, Complex, Complex (0, 0), +); |
|
4569
|
632 } |
|
|
633 |
|
|
634 ComplexNDArray |
|
|
635 ComplexNDArray::prod (int dim) const |
|
|
636 { |
|
|
637 MX_ND_COMPLEX_OP_REDUCTION (*= elem (iter_idx), Complex (1, 0)); |
|
|
638 } |
|
|
639 |
|
|
640 ComplexNDArray |
|
|
641 ComplexNDArray::sumsq (int dim) const |
|
|
642 { |
|
|
643 MX_ND_COMPLEX_OP_REDUCTION |
|
|
644 (+= imag (elem (iter_idx)) |
|
|
645 ? elem (iter_idx) * conj (elem (iter_idx)) |
|
|
646 : std::pow (elem (iter_idx), 2), Complex (0, 0)); |
|
|
647 } |
|
|
648 |
|
|
649 ComplexNDArray |
|
|
650 ComplexNDArray::sum (int dim) const |
|
|
651 { |
|
|
652 MX_ND_COMPLEX_OP_REDUCTION (+= elem (iter_idx), Complex (0, 0)); |
|
|
653 } |
|
|
654 |
|
4915
|
655 ComplexNDArray |
|
5073
|
656 ComplexNDArray::concat (const ComplexNDArray& rb, const Array<int>& ra_idx) |
|
4915
|
657 { |
|
4940
|
658 if (rb.numel () > 0) |
|
5073
|
659 insert (rb, ra_idx); |
|
|
660 return *this; |
|
4915
|
661 } |
|
|
662 |
|
|
663 ComplexNDArray |
|
5073
|
664 ComplexNDArray::concat (const NDArray& rb, const Array<int>& ra_idx) |
|
4758
|
665 { |
|
4915
|
666 ComplexNDArray tmp (rb); |
|
4940
|
667 if (rb.numel () > 0) |
|
5073
|
668 insert (tmp, ra_idx); |
|
|
669 return *this; |
|
4915
|
670 } |
|
|
671 |
|
|
672 ComplexNDArray |
|
5073
|
673 concat (NDArray& ra, ComplexNDArray& rb, const Array<int>& ra_idx) |
|
4915
|
674 { |
|
|
675 ComplexNDArray retval (ra); |
|
4940
|
676 if (rb.numel () > 0) |
|
4915
|
677 retval.insert (rb, ra_idx); |
|
|
678 return retval; |
|
4758
|
679 } |
|
|
680 |
|
4844
|
681 static const Complex Complex_NaN_result (octave_NaN, octave_NaN); |
|
|
682 |
|
|
683 ComplexNDArray |
|
|
684 ComplexNDArray::max (int dim) const |
|
|
685 { |
|
|
686 ArrayN<int> dummy_idx; |
|
|
687 return max (dummy_idx, dim); |
|
|
688 } |
|
|
689 |
|
|
690 ComplexNDArray |
|
|
691 ComplexNDArray::max (ArrayN<int>& idx_arg, int dim) const |
|
|
692 { |
|
|
693 dim_vector dv = dims (); |
|
|
694 dim_vector dr = dims (); |
|
|
695 |
|
|
696 if (dv.numel () == 0 || dim > dv.length () || dim < 0) |
|
|
697 return ComplexNDArray (); |
|
|
698 |
|
|
699 dr(dim) = 1; |
|
|
700 |
|
|
701 ComplexNDArray result (dr); |
|
|
702 idx_arg.resize (dr); |
|
|
703 |
|
|
704 int x_stride = 1; |
|
|
705 int x_len = dv(dim); |
|
|
706 for (int i = 0; i < dim; i++) |
|
|
707 x_stride *= dv(i); |
|
|
708 |
|
|
709 for (int i = 0; i < dr.numel (); i++) |
|
|
710 { |
|
|
711 int x_offset; |
|
|
712 if (x_stride == 1) |
|
|
713 x_offset = i * x_len; |
|
|
714 else |
|
|
715 { |
|
|
716 int x_offset2 = 0; |
|
|
717 x_offset = i; |
|
|
718 while (x_offset >= x_stride) |
|
|
719 { |
|
|
720 x_offset -= x_stride; |
|
|
721 x_offset2++; |
|
|
722 } |
|
|
723 x_offset += x_offset2 * x_stride * x_len; |
|
|
724 } |
|
|
725 |
|
|
726 int idx_j; |
|
|
727 |
|
|
728 Complex tmp_max; |
|
|
729 |
|
|
730 double abs_max = octave_NaN; |
|
|
731 |
|
|
732 for (idx_j = 0; idx_j < x_len; idx_j++) |
|
|
733 { |
|
|
734 tmp_max = elem (idx_j * x_stride + x_offset); |
|
|
735 |
|
|
736 if (! octave_is_NaN_or_NA (tmp_max)) |
|
|
737 { |
|
|
738 abs_max = ::abs(tmp_max); |
|
|
739 break; |
|
|
740 } |
|
|
741 } |
|
|
742 |
|
|
743 for (int j = idx_j+1; j < x_len; j++) |
|
|
744 { |
|
|
745 Complex tmp = elem (j * x_stride + x_offset); |
|
|
746 |
|
|
747 if (octave_is_NaN_or_NA (tmp)) |
|
|
748 continue; |
|
|
749 |
|
|
750 double abs_tmp = ::abs (tmp); |
|
|
751 |
|
|
752 if (abs_tmp > abs_max) |
|
|
753 { |
|
|
754 idx_j = j; |
|
|
755 tmp_max = tmp; |
|
|
756 abs_max = abs_tmp; |
|
|
757 } |
|
|
758 } |
|
|
759 |
|
|
760 if (octave_is_NaN_or_NA (tmp_max)) |
|
|
761 { |
|
|
762 result.elem (i) = Complex_NaN_result; |
|
|
763 idx_arg.elem (i) = 0; |
|
|
764 } |
|
|
765 else |
|
|
766 { |
|
|
767 result.elem (i) = tmp_max; |
|
|
768 idx_arg.elem (i) = idx_j; |
|
|
769 } |
|
|
770 } |
|
|
771 |
|
|
772 return result; |
|
|
773 } |
|
|
774 |
|
|
775 ComplexNDArray |
|
|
776 ComplexNDArray::min (int dim) const |
|
|
777 { |
|
|
778 ArrayN<int> dummy_idx; |
|
|
779 return min (dummy_idx, dim); |
|
|
780 } |
|
|
781 |
|
|
782 ComplexNDArray |
|
|
783 ComplexNDArray::min (ArrayN<int>& idx_arg, int dim) const |
|
|
784 { |
|
|
785 dim_vector dv = dims (); |
|
|
786 dim_vector dr = dims (); |
|
|
787 |
|
|
788 if (dv.numel () == 0 || dim > dv.length () || dim < 0) |
|
|
789 return ComplexNDArray (); |
|
|
790 |
|
|
791 dr(dim) = 1; |
|
|
792 |
|
|
793 ComplexNDArray result (dr); |
|
|
794 idx_arg.resize (dr); |
|
|
795 |
|
|
796 int x_stride = 1; |
|
|
797 int x_len = dv(dim); |
|
|
798 for (int i = 0; i < dim; i++) |
|
|
799 x_stride *= dv(i); |
|
|
800 |
|
|
801 for (int i = 0; i < dr.numel (); i++) |
|
|
802 { |
|
|
803 int x_offset; |
|
|
804 if (x_stride == 1) |
|
|
805 x_offset = i * x_len; |
|
|
806 else |
|
|
807 { |
|
|
808 int x_offset2 = 0; |
|
|
809 x_offset = i; |
|
|
810 while (x_offset >= x_stride) |
|
|
811 { |
|
|
812 x_offset -= x_stride; |
|
|
813 x_offset2++; |
|
|
814 } |
|
|
815 x_offset += x_offset2 * x_stride * x_len; |
|
|
816 } |
|
|
817 |
|
|
818 int idx_j; |
|
|
819 |
|
|
820 Complex tmp_min; |
|
|
821 |
|
|
822 double abs_min = octave_NaN; |
|
|
823 |
|
|
824 for (idx_j = 0; idx_j < x_len; idx_j++) |
|
|
825 { |
|
|
826 tmp_min = elem (idx_j * x_stride + x_offset); |
|
|
827 |
|
|
828 if (! octave_is_NaN_or_NA (tmp_min)) |
|
|
829 { |
|
|
830 abs_min = ::abs(tmp_min); |
|
|
831 break; |
|
|
832 } |
|
|
833 } |
|
|
834 |
|
|
835 for (int j = idx_j+1; j < x_len; j++) |
|
|
836 { |
|
|
837 Complex tmp = elem (j * x_stride + x_offset); |
|
|
838 |
|
|
839 if (octave_is_NaN_or_NA (tmp)) |
|
|
840 continue; |
|
|
841 |
|
|
842 double abs_tmp = ::abs (tmp); |
|
|
843 |
|
|
844 if (abs_tmp < abs_min) |
|
|
845 { |
|
|
846 idx_j = j; |
|
|
847 tmp_min = tmp; |
|
|
848 abs_min = abs_tmp; |
|
|
849 } |
|
|
850 } |
|
|
851 |
|
|
852 if (octave_is_NaN_or_NA (tmp_min)) |
|
|
853 { |
|
|
854 result.elem (i) = Complex_NaN_result; |
|
|
855 idx_arg.elem (i) = 0; |
|
|
856 } |
|
|
857 else |
|
|
858 { |
|
|
859 result.elem (i) = tmp_min; |
|
|
860 idx_arg.elem (i) = idx_j; |
|
|
861 } |
|
|
862 } |
|
|
863 |
|
|
864 return result; |
|
|
865 } |
|
|
866 |
|
4634
|
867 NDArray |
|
4569
|
868 ComplexNDArray::abs (void) const |
|
|
869 { |
|
4634
|
870 NDArray retval (dims ()); |
|
4569
|
871 |
|
4634
|
872 int nel = nelem (); |
|
|
873 |
|
|
874 for (int i = 0; i < nel; i++) |
|
|
875 retval(i) = ::abs (elem (i)); |
|
4569
|
876 |
|
|
877 return retval; |
|
4514
|
878 } |
|
|
879 |
|
4765
|
880 ComplexNDArray& |
|
|
881 ComplexNDArray::insert (const NDArray& a, int r, int c) |
|
|
882 { |
|
|
883 dim_vector a_dv = a.dims (); |
|
|
884 |
|
|
885 int n = a_dv.length (); |
|
|
886 |
|
|
887 if (n == dimensions.length ()) |
|
|
888 { |
|
|
889 Array<int> a_ra_idx (a_dv.length (), 0); |
|
|
890 |
|
|
891 a_ra_idx.elem (0) = r; |
|
|
892 a_ra_idx.elem (1) = c; |
|
|
893 |
|
|
894 for (int i = 0; i < n; i++) |
|
|
895 { |
|
|
896 if (a_ra_idx (i) < 0 || (a_ra_idx (i) + a_dv (i)) > dimensions (i)) |
|
|
897 { |
|
|
898 (*current_liboctave_error_handler) |
|
|
899 ("Array<T>::insert: range error for insert"); |
|
|
900 return *this; |
|
|
901 } |
|
|
902 } |
|
|
903 |
|
|
904 a_ra_idx.elem (0) = 0; |
|
|
905 a_ra_idx.elem (1) = 0; |
|
|
906 |
|
|
907 int n_elt = a.numel (); |
|
|
908 |
|
|
909 // IS make_unique () NECCESSARY HERE?? |
|
|
910 |
|
|
911 for (int i = 0; i < n_elt; i++) |
|
|
912 { |
|
|
913 Array<int> ra_idx = a_ra_idx; |
|
|
914 |
|
|
915 ra_idx.elem (0) = a_ra_idx (0) + r; |
|
|
916 ra_idx.elem (1) = a_ra_idx (1) + c; |
|
|
917 |
|
|
918 elem (ra_idx) = a.elem (a_ra_idx); |
|
|
919 |
|
|
920 increment_index (a_ra_idx, a_dv); |
|
|
921 } |
|
|
922 } |
|
|
923 else |
|
|
924 (*current_liboctave_error_handler) |
|
|
925 ("Array<T>::insert: invalid indexing operation"); |
|
|
926 |
|
|
927 return *this; |
|
|
928 } |
|
|
929 |
|
|
930 ComplexNDArray& |
|
|
931 ComplexNDArray::insert (const ComplexNDArray& a, int r, int c) |
|
|
932 { |
|
|
933 Array<Complex>::insert (a, r, c); |
|
|
934 return *this; |
|
|
935 } |
|
|
936 |
|
4915
|
937 ComplexNDArray& |
|
|
938 ComplexNDArray::insert (const ComplexNDArray& a, const Array<int>& ra_idx) |
|
|
939 { |
|
|
940 Array<Complex>::insert (a, ra_idx); |
|
|
941 return *this; |
|
|
942 } |
|
|
943 |
|
4514
|
944 ComplexMatrix |
|
|
945 ComplexNDArray::matrix_value (void) const |
|
|
946 { |
|
|
947 ComplexMatrix retval; |
|
|
948 |
|
|
949 int nd = ndims (); |
|
|
950 |
|
|
951 switch (nd) |
|
|
952 { |
|
|
953 case 1: |
|
|
954 retval = ComplexMatrix (Array2<Complex> (*this, dimensions(0), 1)); |
|
|
955 break; |
|
|
956 |
|
|
957 case 2: |
|
|
958 retval = ComplexMatrix (Array2<Complex> (*this, dimensions(0), |
|
|
959 dimensions(1))); |
|
|
960 break; |
|
|
961 |
|
|
962 default: |
|
|
963 (*current_liboctave_error_handler) |
|
4770
|
964 ("invalid conversion of ComplexNDArray to ComplexMatrix"); |
|
4514
|
965 break; |
|
|
966 } |
|
|
967 |
|
|
968 return retval; |
|
|
969 } |
|
|
970 |
|
4532
|
971 void |
|
|
972 ComplexNDArray::increment_index (Array<int>& ra_idx, |
|
|
973 const dim_vector& dimensions, |
|
|
974 int start_dimension) |
|
|
975 { |
|
|
976 ::increment_index (ra_idx, dimensions, start_dimension); |
|
|
977 } |
|
|
978 |
|
4556
|
979 int |
|
|
980 ComplexNDArray::compute_index (Array<int>& ra_idx, |
|
|
981 const dim_vector& dimensions) |
|
|
982 { |
|
|
983 return ::compute_index (ra_idx, dimensions); |
|
|
984 } |
|
|
985 |
|
4687
|
986 |
|
|
987 // This contains no information on the array structure !!! |
|
|
988 std::ostream& |
|
|
989 operator << (std::ostream& os, const ComplexNDArray& a) |
|
|
990 { |
|
|
991 int nel = a.nelem (); |
|
|
992 |
|
|
993 for (int i = 0; i < nel; i++) |
|
|
994 { |
|
|
995 os << " "; |
|
|
996 octave_write_complex (os, a.elem (i)); |
|
|
997 os << "\n"; |
|
|
998 } |
|
|
999 return os; |
|
|
1000 } |
|
|
1001 |
|
|
1002 std::istream& |
|
|
1003 operator >> (std::istream& is, ComplexNDArray& a) |
|
|
1004 { |
|
|
1005 int nel = a.nelem (); |
|
|
1006 |
|
|
1007 if (nel < 1 ) |
|
|
1008 is.clear (std::ios::badbit); |
|
|
1009 else |
|
|
1010 { |
|
|
1011 Complex tmp; |
|
|
1012 for (int i = 0; i < nel; i++) |
|
|
1013 { |
|
|
1014 tmp = octave_read_complex (is); |
|
|
1015 if (is) |
|
|
1016 a.elem (i) = tmp; |
|
|
1017 else |
|
|
1018 goto done; |
|
|
1019 } |
|
|
1020 } |
|
|
1021 |
|
|
1022 done: |
|
|
1023 |
|
|
1024 return is; |
|
|
1025 } |
|
|
1026 |
|
4844
|
1027 // XXX FIXME XXX -- it would be nice to share code among the min/max |
|
|
1028 // functions below. |
|
|
1029 |
|
|
1030 #define EMPTY_RETURN_CHECK(T) \ |
|
|
1031 if (nel == 0) \ |
|
|
1032 return T (dv); |
|
|
1033 |
|
|
1034 ComplexNDArray |
|
|
1035 min (const Complex& c, const ComplexNDArray& m) |
|
|
1036 { |
|
|
1037 dim_vector dv = m.dims (); |
|
|
1038 int nel = dv.numel (); |
|
|
1039 |
|
|
1040 EMPTY_RETURN_CHECK (ComplexNDArray); |
|
|
1041 |
|
|
1042 ComplexNDArray result (dv); |
|
|
1043 |
|
|
1044 for (int i = 0; i < nel; i++) |
|
|
1045 { |
|
|
1046 OCTAVE_QUIT; |
|
|
1047 result (i) = xmin (c, m (i)); |
|
|
1048 } |
|
|
1049 |
|
|
1050 return result; |
|
|
1051 } |
|
|
1052 |
|
|
1053 ComplexNDArray |
|
|
1054 min (const ComplexNDArray& m, const Complex& c) |
|
|
1055 { |
|
|
1056 dim_vector dv = m.dims (); |
|
|
1057 int nel = dv.numel (); |
|
|
1058 |
|
|
1059 EMPTY_RETURN_CHECK (ComplexNDArray); |
|
|
1060 |
|
|
1061 ComplexNDArray result (dv); |
|
|
1062 |
|
|
1063 for (int i = 0; i < nel; i++) |
|
|
1064 { |
|
|
1065 OCTAVE_QUIT; |
|
|
1066 result (i) = xmin (c, m (i)); |
|
|
1067 } |
|
|
1068 |
|
|
1069 return result; |
|
|
1070 } |
|
|
1071 |
|
|
1072 ComplexNDArray |
|
|
1073 min (const ComplexNDArray& a, const ComplexNDArray& b) |
|
|
1074 { |
|
|
1075 dim_vector dv = a.dims (); |
|
|
1076 int nel = dv.numel (); |
|
|
1077 |
|
|
1078 if (dv != b.dims ()) |
|
|
1079 { |
|
|
1080 (*current_liboctave_error_handler) |
|
|
1081 ("two-arg min expecting args of same size"); |
|
|
1082 return ComplexNDArray (); |
|
|
1083 } |
|
|
1084 |
|
|
1085 EMPTY_RETURN_CHECK (ComplexNDArray); |
|
|
1086 |
|
|
1087 ComplexNDArray result (dv); |
|
|
1088 |
|
|
1089 for (int i = 0; i < nel; i++) |
|
|
1090 { |
|
|
1091 OCTAVE_QUIT; |
|
|
1092 result (i) = xmin (a (i), b (i)); |
|
|
1093 } |
|
|
1094 |
|
|
1095 return result; |
|
|
1096 } |
|
|
1097 |
|
|
1098 ComplexNDArray |
|
|
1099 max (const Complex& c, const ComplexNDArray& m) |
|
|
1100 { |
|
|
1101 dim_vector dv = m.dims (); |
|
|
1102 int nel = dv.numel (); |
|
|
1103 |
|
|
1104 EMPTY_RETURN_CHECK (ComplexNDArray); |
|
|
1105 |
|
|
1106 ComplexNDArray result (dv); |
|
|
1107 |
|
|
1108 for (int i = 0; i < nel; i++) |
|
|
1109 { |
|
|
1110 OCTAVE_QUIT; |
|
|
1111 result (i) = xmax (c, m (i)); |
|
|
1112 } |
|
|
1113 |
|
|
1114 return result; |
|
|
1115 } |
|
|
1116 |
|
|
1117 ComplexNDArray |
|
|
1118 max (const ComplexNDArray& m, const Complex& c) |
|
|
1119 { |
|
|
1120 dim_vector dv = m.dims (); |
|
|
1121 int nel = dv.numel (); |
|
|
1122 |
|
|
1123 EMPTY_RETURN_CHECK (ComplexNDArray); |
|
|
1124 |
|
|
1125 ComplexNDArray result (dv); |
|
|
1126 |
|
|
1127 for (int i = 0; i < nel; i++) |
|
|
1128 { |
|
|
1129 OCTAVE_QUIT; |
|
|
1130 result (i) = xmax (c, m (i)); |
|
|
1131 } |
|
|
1132 |
|
|
1133 return result; |
|
|
1134 } |
|
|
1135 |
|
|
1136 ComplexNDArray |
|
|
1137 max (const ComplexNDArray& a, const ComplexNDArray& b) |
|
|
1138 { |
|
|
1139 dim_vector dv = a.dims (); |
|
|
1140 int nel = dv.numel (); |
|
|
1141 |
|
|
1142 if (dv != b.dims ()) |
|
|
1143 { |
|
|
1144 (*current_liboctave_error_handler) |
|
|
1145 ("two-arg max expecting args of same size"); |
|
|
1146 return ComplexNDArray (); |
|
|
1147 } |
|
|
1148 |
|
|
1149 EMPTY_RETURN_CHECK (ComplexNDArray); |
|
|
1150 |
|
|
1151 ComplexNDArray result (dv); |
|
|
1152 |
|
|
1153 for (int i = 0; i < nel; i++) |
|
|
1154 { |
|
|
1155 OCTAVE_QUIT; |
|
|
1156 result (i) = xmax (a (i), b (i)); |
|
|
1157 } |
|
|
1158 |
|
|
1159 return result; |
|
|
1160 } |
|
|
1161 |
|
4543
|
1162 NDS_CMP_OPS(ComplexNDArray, real, Complex, real) |
|
|
1163 NDS_BOOL_OPS(ComplexNDArray, Complex, 0.0) |
|
|
1164 |
|
|
1165 SND_CMP_OPS(Complex, real, ComplexNDArray, real) |
|
|
1166 SND_BOOL_OPS(Complex, ComplexNDArray, 0.0) |
|
|
1167 |
|
|
1168 NDND_CMP_OPS(ComplexNDArray, real, ComplexNDArray, real) |
|
|
1169 NDND_BOOL_OPS(ComplexNDArray, ComplexNDArray, 0.0) |
|
|
1170 |
|
4514
|
1171 /* |
|
|
1172 ;;; Local Variables: *** |
|
|
1173 ;;; mode: C++ *** |
|
|
1174 ;;; End: *** |
|
|
1175 */ |
