1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1996-2021 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29
30 #include <complex>
31 #include <istream>
32 #include <ostream>
33
34 #include "Array-util.h"
35 #include "f77-fcn.h"
36 #include "fCNDArray.h"
37 #include "lo-ieee.h"
38 #include "lo-mappers.h"
39 #include "mx-base.h"
40 #include "mx-op-defs.h"
41 #include "mx-fcnda-fs.h"
42 #include "oct-fftw.h"
43 #include "oct-locbuf.h"
44
45 #include "bsxfun-defs.cc"
46
FloatComplexNDArray(const charNDArray & a)47 FloatComplexNDArray::FloatComplexNDArray (const charNDArray& a)
48 : MArray<FloatComplex> (a.dims ())
49 {
50 octave_idx_type n = a.numel ();
51 for (octave_idx_type i = 0; i < n; i++)
52 xelem (i) = static_cast<unsigned char> (a(i));
53 }
54
55 #if defined (HAVE_FFTW)
56
57 FloatComplexNDArray
fourier(int dim) const58 FloatComplexNDArray::fourier (int dim) const
59 {
60 dim_vector dv = dims ();
61
62 if (dim > dv.ndims () || dim < 0)
63 return FloatComplexNDArray ();
64
65 octave_idx_type stride = 1;
66 octave_idx_type n = dv(dim);
67
68 for (int i = 0; i < dim; i++)
69 stride *= dv(i);
70
71 octave_idx_type howmany = numel () / dv(dim);
72 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
73 octave_idx_type nloop = (stride == 1 ? 1 : numel () / dv(dim) / stride);
74 octave_idx_type dist = (stride == 1 ? n : 1);
75
76 const FloatComplex *in (fortran_vec ());
77 FloatComplexNDArray retval (dv);
78 FloatComplex *out (retval.fortran_vec ());
79
80 // Need to be careful here about the distance between fft's
81 for (octave_idx_type k = 0; k < nloop; k++)
82 octave::fftw::fft (in + k * stride * n, out + k * stride * n,
83 n, howmany, stride, dist);
84
85 return retval;
86 }
87
88 FloatComplexNDArray
ifourier(int dim) const89 FloatComplexNDArray::ifourier (int dim) const
90 {
91 dim_vector dv = dims ();
92
93 if (dim > dv.ndims () || dim < 0)
94 return FloatComplexNDArray ();
95
96 octave_idx_type stride = 1;
97 octave_idx_type n = dv(dim);
98
99 for (int i = 0; i < dim; i++)
100 stride *= dv(i);
101
102 octave_idx_type howmany = numel () / dv(dim);
103 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
104 octave_idx_type nloop = (stride == 1 ? 1 : numel () / dv(dim) / stride);
105 octave_idx_type dist = (stride == 1 ? n : 1);
106
107 const FloatComplex *in (fortran_vec ());
108 FloatComplexNDArray retval (dv);
109 FloatComplex *out (retval.fortran_vec ());
110
111 // Need to be careful here about the distance between fft's
112 for (octave_idx_type k = 0; k < nloop; k++)
113 octave::fftw::ifft (in + k * stride * n, out + k * stride * n,
114 n, howmany, stride, dist);
115
116 return retval;
117 }
118
119 FloatComplexNDArray
fourier2d(void) const120 FloatComplexNDArray::fourier2d (void) const
121 {
122 dim_vector dv = dims ();
123 if (dv.ndims () < 2)
124 return FloatComplexNDArray ();
125
126 dim_vector dv2 (dv(0), dv(1));
127 const FloatComplex *in = fortran_vec ();
128 FloatComplexNDArray retval (dv);
129 FloatComplex *out = retval.fortran_vec ();
130 octave_idx_type howmany = numel () / dv(0) / dv(1);
131 octave_idx_type dist = dv(0) * dv(1);
132
133 for (octave_idx_type i=0; i < howmany; i++)
134 octave::fftw::fftNd (in + i*dist, out + i*dist, 2, dv2);
135
136 return retval;
137 }
138
139 FloatComplexNDArray
ifourier2d(void) const140 FloatComplexNDArray::ifourier2d (void) const
141 {
142 dim_vector dv = dims ();
143 if (dv.ndims () < 2)
144 return FloatComplexNDArray ();
145
146 dim_vector dv2 (dv(0), dv(1));
147 const FloatComplex *in = fortran_vec ();
148 FloatComplexNDArray retval (dv);
149 FloatComplex *out = retval.fortran_vec ();
150 octave_idx_type howmany = numel () / dv(0) / dv(1);
151 octave_idx_type dist = dv(0) * dv(1);
152
153 for (octave_idx_type i=0; i < howmany; i++)
154 octave::fftw::ifftNd (in + i*dist, out + i*dist, 2, dv2);
155
156 return retval;
157 }
158
159 FloatComplexNDArray
fourierNd(void) const160 FloatComplexNDArray::fourierNd (void) const
161 {
162 dim_vector dv = dims ();
163 int rank = dv.ndims ();
164
165 const FloatComplex *in (fortran_vec ());
166 FloatComplexNDArray retval (dv);
167 FloatComplex *out (retval.fortran_vec ());
168
169 octave::fftw::fftNd (in, out, rank, dv);
170
171 return retval;
172 }
173
174 FloatComplexNDArray
ifourierNd(void) const175 FloatComplexNDArray::ifourierNd (void) const
176 {
177 dim_vector dv = dims ();
178 int rank = dv.ndims ();
179
180 const FloatComplex *in (fortran_vec ());
181 FloatComplexNDArray retval (dv);
182 FloatComplex *out (retval.fortran_vec ());
183
184 octave::fftw::ifftNd (in, out, rank, dv);
185
186 return retval;
187 }
188
189 #else
190
191 FloatComplexNDArray
fourier(int dim) const192 FloatComplexNDArray::fourier (int dim) const
193 {
194 octave_unused_parameter (dim);
195
196 (*current_liboctave_error_handler)
197 ("support for FFTW was unavailable or disabled when liboctave was built");
198
199 return FloatComplexNDArray ();
200 }
201
202 FloatComplexNDArray
ifourier(int dim) const203 FloatComplexNDArray::ifourier (int dim) const
204 {
205 octave_unused_parameter (dim);
206
207 (*current_liboctave_error_handler)
208 ("support for FFTW was unavailable or disabled when liboctave was built");
209
210 return FloatComplexNDArray ();
211 }
212
213 FloatComplexNDArray
fourier2d(void) const214 FloatComplexNDArray::fourier2d (void) const
215 {
216 (*current_liboctave_error_handler)
217 ("support for FFTW was unavailable or disabled when liboctave was built");
218
219 return FloatComplexNDArray ();
220 }
221
222 FloatComplexNDArray
ifourier2d(void) const223 FloatComplexNDArray::ifourier2d (void) const
224 {
225 (*current_liboctave_error_handler)
226 ("support for FFTW was unavailable or disabled when liboctave was built");
227
228 return FloatComplexNDArray ();
229 }
230
231 FloatComplexNDArray
fourierNd(void) const232 FloatComplexNDArray::fourierNd (void) const
233 {
234 (*current_liboctave_error_handler)
235 ("support for FFTW was unavailable or disabled when liboctave was built");
236
237 return FloatComplexNDArray ();
238 }
239
240 FloatComplexNDArray
ifourierNd(void) const241 FloatComplexNDArray::ifourierNd (void) const
242 {
243 (*current_liboctave_error_handler)
244 ("support for FFTW was unavailable or disabled when liboctave was built");
245
246 return FloatComplexNDArray ();
247 }
248
249 #endif
250
251 // unary operations
252
253 boolNDArray
operator !(void) const254 FloatComplexNDArray::operator ! (void) const
255 {
256 if (any_element_is_nan ())
257 octave::err_nan_to_logical_conversion ();
258
259 return do_mx_unary_op<bool, FloatComplex> (*this, mx_inline_not);
260 }
261
262 // FIXME: this is not quite the right thing.
263
264 bool
any_element_is_nan(void) const265 FloatComplexNDArray::any_element_is_nan (void) const
266 {
267 return do_mx_check<FloatComplex> (*this, mx_inline_any_nan);
268 }
269
270 bool
any_element_is_inf_or_nan(void) const271 FloatComplexNDArray::any_element_is_inf_or_nan (void) const
272 {
273 return ! do_mx_check<FloatComplex> (*this, mx_inline_all_finite);
274 }
275
276 // Return true if no elements have imaginary components.
277
278 bool
all_elements_are_real(void) const279 FloatComplexNDArray::all_elements_are_real (void) const
280 {
281 return do_mx_check<FloatComplex> (*this, mx_inline_all_real);
282 }
283
284 // Return nonzero if any element of CM has a non-integer real or
285 // imaginary part. Also extract the largest and smallest (real or
286 // imaginary) values and return them in MAX_VAL and MIN_VAL.
287
288 bool
all_integers(float & max_val,float & min_val) const289 FloatComplexNDArray::all_integers (float& max_val, float& min_val) const
290 {
291 octave_idx_type nel = numel ();
292
293 if (nel > 0)
294 {
295 FloatComplex val = elem (0);
296
297 float r_val = val.real ();
298 float i_val = val.imag ();
299
300 max_val = r_val;
301 min_val = r_val;
302
303 if (i_val > max_val)
304 max_val = i_val;
305
306 if (i_val < max_val)
307 min_val = i_val;
308 }
309 else
310 return false;
311
312 for (octave_idx_type i = 0; i < nel; i++)
313 {
314 FloatComplex val = elem (i);
315
316 float r_val = val.real ();
317 float i_val = val.imag ();
318
319 if (r_val > max_val)
320 max_val = r_val;
321
322 if (i_val > max_val)
323 max_val = i_val;
324
325 if (r_val < min_val)
326 min_val = r_val;
327
328 if (i_val < min_val)
329 min_val = i_val;
330
331 if (octave::math::x_nint (r_val) != r_val
332 || octave::math::x_nint (i_val) != i_val)
333 return false;
334 }
335
336 return true;
337 }
338
339 bool
too_large_for_float(void) const340 FloatComplexNDArray::too_large_for_float (void) const
341 {
342 return false;
343 }
344
345 boolNDArray
all(int dim) const346 FloatComplexNDArray::all (int dim) const
347 {
348 return do_mx_red_op<bool, FloatComplex> (*this, dim, mx_inline_all);
349 }
350
351 boolNDArray
any(int dim) const352 FloatComplexNDArray::any (int dim) const
353 {
354 return do_mx_red_op<bool, FloatComplex> (*this, dim, mx_inline_any);
355 }
356
357 FloatComplexNDArray
cumprod(int dim) const358 FloatComplexNDArray::cumprod (int dim) const
359 {
360 return do_mx_cum_op<FloatComplex, FloatComplex> (*this, dim,
361 mx_inline_cumprod);
362 }
363
364 FloatComplexNDArray
cumsum(int dim) const365 FloatComplexNDArray::cumsum (int dim) const
366 {
367 return do_mx_cum_op<FloatComplex, FloatComplex> (*this, dim,
368 mx_inline_cumsum);
369 }
370
371 FloatComplexNDArray
prod(int dim) const372 FloatComplexNDArray::prod (int dim) const
373 {
374 return do_mx_red_op<FloatComplex, FloatComplex> (*this, dim, mx_inline_prod);
375 }
376
377 ComplexNDArray
dprod(int dim) const378 FloatComplexNDArray::dprod (int dim) const
379 {
380 return do_mx_red_op<Complex, FloatComplex> (*this, dim, mx_inline_dprod);
381 }
382
383 FloatComplexNDArray
sum(int dim) const384 FloatComplexNDArray::sum (int dim) const
385 {
386 return do_mx_red_op<FloatComplex, FloatComplex> (*this, dim, mx_inline_sum);
387 }
388
389 ComplexNDArray
dsum(int dim) const390 FloatComplexNDArray::dsum (int dim) const
391 {
392 return do_mx_red_op<Complex, FloatComplex> (*this, dim, mx_inline_dsum);
393 }
394
395 FloatComplexNDArray
sumsq(int dim) const396 FloatComplexNDArray::sumsq (int dim) const
397 {
398 return do_mx_red_op<float, FloatComplex> (*this, dim, mx_inline_sumsq);
399 }
400
401 FloatComplexNDArray
diff(octave_idx_type order,int dim) const402 FloatComplexNDArray::diff (octave_idx_type order, int dim) const
403 {
404 return do_mx_diff_op<FloatComplex> (*this, dim, order, mx_inline_diff);
405 }
406
407 FloatComplexNDArray
concat(const FloatComplexNDArray & rb,const Array<octave_idx_type> & ra_idx)408 FloatComplexNDArray::concat (const FloatComplexNDArray& rb,
409 const Array<octave_idx_type>& ra_idx)
410 {
411 if (rb.numel () > 0)
412 insert (rb, ra_idx);
413 return *this;
414 }
415
416 FloatComplexNDArray
concat(const FloatNDArray & rb,const Array<octave_idx_type> & ra_idx)417 FloatComplexNDArray::concat (const FloatNDArray& rb,
418 const Array<octave_idx_type>& ra_idx)
419 {
420 FloatComplexNDArray tmp (rb);
421 if (rb.numel () > 0)
422 insert (tmp, ra_idx);
423 return *this;
424 }
425
426 FloatComplexNDArray
concat(NDArray & ra,FloatComplexNDArray & rb,const Array<octave_idx_type> & ra_idx)427 concat (NDArray& ra, FloatComplexNDArray& rb,
428 const Array<octave_idx_type>& ra_idx)
429 {
430 FloatComplexNDArray retval (ra);
431 if (rb.numel () > 0)
432 retval.insert (rb, ra_idx);
433 return retval;
434 }
435
436 static const FloatComplex FloatComplex_NaN_result (octave::numeric_limits<float>::NaN (),
437 octave::numeric_limits<float>::NaN ());
438
439 FloatComplexNDArray
max(int dim) const440 FloatComplexNDArray::max (int dim) const
441 {
442 return do_mx_minmax_op<FloatComplex> (*this, dim, mx_inline_max);
443 }
444
445 FloatComplexNDArray
max(Array<octave_idx_type> & idx_arg,int dim) const446 FloatComplexNDArray::max (Array<octave_idx_type>& idx_arg, int dim) const
447 {
448 return do_mx_minmax_op<FloatComplex> (*this, idx_arg, dim, mx_inline_max);
449 }
450
451 FloatComplexNDArray
min(int dim) const452 FloatComplexNDArray::min (int dim) const
453 {
454 return do_mx_minmax_op<FloatComplex> (*this, dim, mx_inline_min);
455 }
456
457 FloatComplexNDArray
min(Array<octave_idx_type> & idx_arg,int dim) const458 FloatComplexNDArray::min (Array<octave_idx_type>& idx_arg, int dim) const
459 {
460 return do_mx_minmax_op<FloatComplex> (*this, idx_arg, dim, mx_inline_min);
461 }
462
463 FloatComplexNDArray
cummax(int dim) const464 FloatComplexNDArray::cummax (int dim) const
465 {
466 return do_mx_cumminmax_op<FloatComplex> (*this, dim, mx_inline_cummax);
467 }
468
469 FloatComplexNDArray
cummax(Array<octave_idx_type> & idx_arg,int dim) const470 FloatComplexNDArray::cummax (Array<octave_idx_type>& idx_arg, int dim) const
471 {
472 return do_mx_cumminmax_op<FloatComplex> (*this, idx_arg, dim,
473 mx_inline_cummax);
474 }
475
476 FloatComplexNDArray
cummin(int dim) const477 FloatComplexNDArray::cummin (int dim) const
478 {
479 return do_mx_cumminmax_op<FloatComplex> (*this, dim, mx_inline_cummin);
480 }
481
482 FloatComplexNDArray
cummin(Array<octave_idx_type> & idx_arg,int dim) const483 FloatComplexNDArray::cummin (Array<octave_idx_type>& idx_arg, int dim) const
484 {
485 return do_mx_cumminmax_op<FloatComplex> (*this, idx_arg, dim,
486 mx_inline_cummin);
487 }
488
489 FloatNDArray
abs(void) const490 FloatComplexNDArray::abs (void) const
491 {
492 return do_mx_unary_map<float, FloatComplex, std::abs> (*this);
493 }
494
495 boolNDArray
isnan(void) const496 FloatComplexNDArray::isnan (void) const
497 {
498 return do_mx_unary_map<bool, FloatComplex, octave::math::isnan> (*this);
499 }
500
501 boolNDArray
isinf(void) const502 FloatComplexNDArray::isinf (void) const
503 {
504 return do_mx_unary_map<bool, FloatComplex, octave::math::isinf> (*this);
505 }
506
507 boolNDArray
isfinite(void) const508 FloatComplexNDArray::isfinite (void) const
509 {
510 return do_mx_unary_map<bool, FloatComplex, octave::math::isfinite> (*this);
511 }
512
513 FloatComplexNDArray
conj(const FloatComplexNDArray & a)514 conj (const FloatComplexNDArray& a)
515 {
516 return do_mx_unary_map<FloatComplex, FloatComplex, std::conj<float>> (a);
517 }
518
519 FloatComplexNDArray&
insert(const NDArray & a,octave_idx_type r,octave_idx_type c)520 FloatComplexNDArray::insert (const NDArray& a,
521 octave_idx_type r, octave_idx_type c)
522 {
523 dim_vector a_dv = a.dims ();
524
525 int n = a_dv.ndims ();
526
527 if (n == dimensions.ndims ())
528 {
529 Array<octave_idx_type> a_ra_idx (dim_vector (a_dv.ndims (), 1), 0);
530
531 a_ra_idx.elem (0) = r;
532 a_ra_idx.elem (1) = c;
533
534 for (int i = 0; i < n; i++)
535 {
536 if (a_ra_idx(i) < 0 || (a_ra_idx(i) + a_dv(i)) > dimensions(i))
537 (*current_liboctave_error_handler)
538 ("Array<T>::insert: range error for insert");
539 }
540
541 a_ra_idx.elem (0) = 0;
542 a_ra_idx.elem (1) = 0;
543
544 octave_idx_type n_elt = a.numel ();
545
546 // IS make_unique () NECESSARY HERE?
547
548 for (octave_idx_type i = 0; i < n_elt; i++)
549 {
550 Array<octave_idx_type> ra_idx = a_ra_idx;
551
552 ra_idx.elem (0) = a_ra_idx(0) + r;
553 ra_idx.elem (1) = a_ra_idx(1) + c;
554
555 elem (ra_idx) = a.elem (a_ra_idx);
556
557 increment_index (a_ra_idx, a_dv);
558 }
559 }
560 else
561 (*current_liboctave_error_handler)
562 ("Array<T>::insert: invalid indexing operation");
563
564 return *this;
565 }
566
567 FloatComplexNDArray&
insert(const FloatComplexNDArray & a,octave_idx_type r,octave_idx_type c)568 FloatComplexNDArray::insert (const FloatComplexNDArray& a,
569 octave_idx_type r, octave_idx_type c)
570 {
571 Array<FloatComplex>::insert (a, r, c);
572 return *this;
573 }
574
575 FloatComplexNDArray&
insert(const FloatComplexNDArray & a,const Array<octave_idx_type> & ra_idx)576 FloatComplexNDArray::insert (const FloatComplexNDArray& a,
577 const Array<octave_idx_type>& ra_idx)
578 {
579 Array<FloatComplex>::insert (a, ra_idx);
580 return *this;
581 }
582
583 void
increment_index(Array<octave_idx_type> & ra_idx,const dim_vector & dimensions,int start_dimension)584 FloatComplexNDArray::increment_index (Array<octave_idx_type>& ra_idx,
585 const dim_vector& dimensions,
586 int start_dimension)
587 {
588 ::increment_index (ra_idx, dimensions, start_dimension);
589 }
590
591 octave_idx_type
compute_index(Array<octave_idx_type> & ra_idx,const dim_vector & dimensions)592 FloatComplexNDArray::compute_index (Array<octave_idx_type>& ra_idx,
593 const dim_vector& dimensions)
594 {
595 return ::compute_index (ra_idx, dimensions);
596 }
597
598 FloatComplexNDArray
diag(octave_idx_type k) const599 FloatComplexNDArray::diag (octave_idx_type k) const
600 {
601 return MArray<FloatComplex>::diag (k);
602 }
603
604 FloatComplexNDArray
diag(octave_idx_type m,octave_idx_type n) const605 FloatComplexNDArray::diag (octave_idx_type m, octave_idx_type n) const
606 {
607 return MArray<FloatComplex>::diag (m, n);
608 }
609
610 // This contains no information on the array structure !!!
611 std::ostream&
operator <<(std::ostream & os,const FloatComplexNDArray & a)612 operator << (std::ostream& os, const FloatComplexNDArray& a)
613 {
614 octave_idx_type nel = a.numel ();
615
616 for (octave_idx_type i = 0; i < nel; i++)
617 {
618 os << ' ';
619 octave_write_complex (os, a.elem (i));
620 os << "\n";
621 }
622 return os;
623 }
624
625 std::istream&
operator >>(std::istream & is,FloatComplexNDArray & a)626 operator >> (std::istream& is, FloatComplexNDArray& a)
627 {
628 octave_idx_type nel = a.numel ();
629
630 if (nel > 0)
631 {
632 FloatComplex tmp;
633 for (octave_idx_type i = 0; i < nel; i++)
634 {
635 tmp = octave_read_value<FloatComplex> (is);
636 if (is)
637 a.elem (i) = tmp;
638 else
639 return is;
640 }
641 }
642
643 return is;
644 }
645
MINMAX_FCNS(FloatComplexNDArray,FloatComplex)646 MINMAX_FCNS (FloatComplexNDArray, FloatComplex)
647
648 NDS_CMP_OPS (FloatComplexNDArray, FloatComplex)
649 NDS_BOOL_OPS (FloatComplexNDArray, FloatComplex)
650
651 SND_CMP_OPS (FloatComplex, FloatComplexNDArray)
652 SND_BOOL_OPS (FloatComplex, FloatComplexNDArray)
653
654 NDND_CMP_OPS (FloatComplexNDArray, FloatComplexNDArray)
655 NDND_BOOL_OPS (FloatComplexNDArray, FloatComplexNDArray)
656
657 FloatComplexNDArray& operator *= (FloatComplexNDArray& a, float s)
658 {
659 if (a.is_shared ())
660 a = a * s;
661 else
662 do_ms_inplace_op<FloatComplex, float> (a, s, mx_inline_mul2);
663 return a;
664 }
665
operator /=(FloatComplexNDArray & a,float s)666 FloatComplexNDArray& operator /= (FloatComplexNDArray& a, float s)
667 {
668 if (a.is_shared ())
669 a = a / s;
670 else
671 do_ms_inplace_op<FloatComplex, float> (a, s, mx_inline_div2);
672 return a;
673 }
674
675 BSXFUN_STDOP_DEFS_MXLOOP (FloatComplexNDArray)
676 BSXFUN_STDREL_DEFS_MXLOOP (FloatComplexNDArray)
677
678 BSXFUN_OP_DEF_MXLOOP (pow, FloatComplexNDArray, mx_inline_pow)
679