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 "CNDArray.h"
36 #include "f77-fcn.h"
37 #include "lo-ieee.h"
38 #include "lo-mappers.h"
39 #include "mx-base.h"
40 #include "mx-cnda-s.h"
41 #include "mx-op-defs.h"
42 #include "oct-fftw.h"
43 #include "oct-locbuf.h"
44
45 #include "bsxfun-defs.cc"
46
ComplexNDArray(const charNDArray & a)47 ComplexNDArray::ComplexNDArray (const charNDArray& a)
48 : MArray<Complex> (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 ComplexNDArray
fourier(int dim) const58 ComplexNDArray::fourier (int dim) const
59 {
60 dim_vector dv = dims ();
61
62 if (dim > dv.ndims () || dim < 0)
63 return ComplexNDArray ();
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 Complex *in (fortran_vec ());
77 ComplexNDArray retval (dv);
78 Complex *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 ComplexNDArray
ifourier(int dim) const89 ComplexNDArray::ifourier (int dim) const
90 {
91 dim_vector dv = dims ();
92
93 if (dim > dv.ndims () || dim < 0)
94 return ComplexNDArray ();
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 Complex *in (fortran_vec ());
108 ComplexNDArray retval (dv);
109 Complex *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 ComplexNDArray
fourier2d(void) const120 ComplexNDArray::fourier2d (void) const
121 {
122 dim_vector dv = dims ();
123 if (dv.ndims () < 2)
124 return ComplexNDArray ();
125
126 dim_vector dv2 (dv(0), dv(1));
127 const Complex *in = fortran_vec ();
128 ComplexNDArray retval (dv);
129 Complex *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 ComplexNDArray
ifourier2d(void) const140 ComplexNDArray::ifourier2d (void) const
141 {
142 dim_vector dv = dims ();
143 if (dv.ndims () < 2)
144 return ComplexNDArray ();
145
146 dim_vector dv2 (dv(0), dv(1));
147 const Complex *in = fortran_vec ();
148 ComplexNDArray retval (dv);
149 Complex *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 ComplexNDArray
fourierNd(void) const160 ComplexNDArray::fourierNd (void) const
161 {
162 dim_vector dv = dims ();
163 int rank = dv.ndims ();
164
165 const Complex *in (fortran_vec ());
166 ComplexNDArray retval (dv);
167 Complex *out (retval.fortran_vec ());
168
169 octave::fftw::fftNd (in, out, rank, dv);
170
171 return retval;
172 }
173
174 ComplexNDArray
ifourierNd(void) const175 ComplexNDArray::ifourierNd (void) const
176 {
177 dim_vector dv = dims ();
178 int rank = dv.ndims ();
179
180 const Complex *in (fortran_vec ());
181 ComplexNDArray retval (dv);
182 Complex *out (retval.fortran_vec ());
183
184 octave::fftw::ifftNd (in, out, rank, dv);
185
186 return retval;
187 }
188
189 #else
190
191 ComplexNDArray
fourier(int dim) const192 ComplexNDArray::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 ComplexNDArray ();
200 }
201
202 ComplexNDArray
ifourier(int dim) const203 ComplexNDArray::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 ComplexNDArray ();
211 }
212
213 ComplexNDArray
fourier2d(void) const214 ComplexNDArray::fourier2d (void) const
215 {
216 (*current_liboctave_error_handler)
217 ("support for FFTW was unavailable or disabled when liboctave was built");
218
219 return ComplexNDArray ();
220 }
221
222 ComplexNDArray
ifourier2d(void) const223 ComplexNDArray::ifourier2d (void) const
224 {
225 (*current_liboctave_error_handler)
226 ("support for FFTW was unavailable or disabled when liboctave was built");
227
228 return ComplexNDArray ();
229 }
230
231 ComplexNDArray
fourierNd(void) const232 ComplexNDArray::fourierNd (void) const
233 {
234 (*current_liboctave_error_handler)
235 ("support for FFTW was unavailable or disabled when liboctave was built");
236
237 return ComplexNDArray ();
238 }
239
240 ComplexNDArray
ifourierNd(void) const241 ComplexNDArray::ifourierNd (void) const
242 {
243 (*current_liboctave_error_handler)
244 ("support for FFTW was unavailable or disabled when liboctave was built");
245
246 return ComplexNDArray ();
247 }
248
249 #endif
250
251 // unary operations
252
253 boolNDArray
operator !(void) const254 ComplexNDArray::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, Complex> (*this, mx_inline_not);
260 }
261
262 // FIXME: this is not quite the right thing.
263
264 bool
any_element_is_nan(void) const265 ComplexNDArray::any_element_is_nan (void) const
266 {
267 return do_mx_check<Complex> (*this, mx_inline_any_nan);
268 }
269
270 bool
any_element_is_inf_or_nan(void) const271 ComplexNDArray::any_element_is_inf_or_nan (void) const
272 {
273 return ! do_mx_check<Complex> (*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 ComplexNDArray::all_elements_are_real (void) const
280 {
281 return do_mx_check<Complex> (*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(double & max_val,double & min_val) const289 ComplexNDArray::all_integers (double& max_val, double& min_val) const
290 {
291 octave_idx_type nel = numel ();
292
293 if (nel > 0)
294 {
295 Complex val = elem (0);
296
297 double r_val = val.real ();
298 double 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 Complex val = elem (i);
315
316 double r_val = val.real ();
317 double 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 ComplexNDArray::too_large_for_float (void) const
341 {
342 return test_any (xtoo_large_for_float);
343 }
344
345 boolNDArray
all(int dim) const346 ComplexNDArray::all (int dim) const
347 {
348 return do_mx_red_op<bool, Complex> (*this, dim, mx_inline_all);
349 }
350
351 boolNDArray
any(int dim) const352 ComplexNDArray::any (int dim) const
353 {
354 return do_mx_red_op<bool, Complex> (*this, dim, mx_inline_any);
355 }
356
357 ComplexNDArray
cumprod(int dim) const358 ComplexNDArray::cumprod (int dim) const
359 {
360 return do_mx_cum_op<Complex, Complex> (*this, dim, mx_inline_cumprod);
361 }
362
363 ComplexNDArray
cumsum(int dim) const364 ComplexNDArray::cumsum (int dim) const
365 {
366 return do_mx_cum_op<Complex, Complex> (*this, dim, mx_inline_cumsum);
367 }
368
369 ComplexNDArray
prod(int dim) const370 ComplexNDArray::prod (int dim) const
371 {
372 return do_mx_red_op<Complex, Complex> (*this, dim, mx_inline_prod);
373 }
374
375 ComplexNDArray
sum(int dim) const376 ComplexNDArray::sum (int dim) const
377 {
378 return do_mx_red_op<Complex, Complex> (*this, dim, mx_inline_sum);
379 }
380
381 ComplexNDArray
xsum(int dim) const382 ComplexNDArray::xsum (int dim) const
383 {
384 return do_mx_red_op<Complex, Complex> (*this, dim, mx_inline_xsum);
385 }
386
387 ComplexNDArray
sumsq(int dim) const388 ComplexNDArray::sumsq (int dim) const
389 {
390 return do_mx_red_op<double, Complex> (*this, dim, mx_inline_sumsq);
391 }
392
393 ComplexNDArray
diff(octave_idx_type order,int dim) const394 ComplexNDArray::diff (octave_idx_type order, int dim) const
395 {
396 return do_mx_diff_op<Complex> (*this, dim, order, mx_inline_diff);
397 }
398
399 ComplexNDArray
concat(const ComplexNDArray & rb,const Array<octave_idx_type> & ra_idx)400 ComplexNDArray::concat (const ComplexNDArray& rb,
401 const Array<octave_idx_type>& ra_idx)
402 {
403 if (rb.numel () > 0)
404 insert (rb, ra_idx);
405 return *this;
406 }
407
408 ComplexNDArray
concat(const NDArray & rb,const Array<octave_idx_type> & ra_idx)409 ComplexNDArray::concat (const NDArray& rb, const Array<octave_idx_type>& ra_idx)
410 {
411 ComplexNDArray tmp (rb);
412 if (rb.numel () > 0)
413 insert (tmp, ra_idx);
414 return *this;
415 }
416
417 ComplexNDArray
concat(NDArray & ra,ComplexNDArray & rb,const Array<octave_idx_type> & ra_idx)418 concat (NDArray& ra, ComplexNDArray& rb, const Array<octave_idx_type>& ra_idx)
419 {
420 ComplexNDArray retval (ra);
421 if (rb.numel () > 0)
422 retval.insert (rb, ra_idx);
423 return retval;
424 }
425
426 static const Complex Complex_NaN_result (octave::numeric_limits<double>::NaN (),
427 octave::numeric_limits<double>::NaN ());
428
429 ComplexNDArray
max(int dim) const430 ComplexNDArray::max (int dim) const
431 {
432 return do_mx_minmax_op<Complex> (*this, dim, mx_inline_max);
433 }
434
435 ComplexNDArray
max(Array<octave_idx_type> & idx_arg,int dim) const436 ComplexNDArray::max (Array<octave_idx_type>& idx_arg, int dim) const
437 {
438 return do_mx_minmax_op<Complex> (*this, idx_arg, dim, mx_inline_max);
439 }
440
441 ComplexNDArray
min(int dim) const442 ComplexNDArray::min (int dim) const
443 {
444 return do_mx_minmax_op<Complex> (*this, dim, mx_inline_min);
445 }
446
447 ComplexNDArray
min(Array<octave_idx_type> & idx_arg,int dim) const448 ComplexNDArray::min (Array<octave_idx_type>& idx_arg, int dim) const
449 {
450 return do_mx_minmax_op<Complex> (*this, idx_arg, dim, mx_inline_min);
451 }
452
453 ComplexNDArray
cummax(int dim) const454 ComplexNDArray::cummax (int dim) const
455 {
456 return do_mx_cumminmax_op<Complex> (*this, dim, mx_inline_cummax);
457 }
458
459 ComplexNDArray
cummax(Array<octave_idx_type> & idx_arg,int dim) const460 ComplexNDArray::cummax (Array<octave_idx_type>& idx_arg, int dim) const
461 {
462 return do_mx_cumminmax_op<Complex> (*this, idx_arg, dim, mx_inline_cummax);
463 }
464
465 ComplexNDArray
cummin(int dim) const466 ComplexNDArray::cummin (int dim) const
467 {
468 return do_mx_cumminmax_op<Complex> (*this, dim, mx_inline_cummin);
469 }
470
471 ComplexNDArray
cummin(Array<octave_idx_type> & idx_arg,int dim) const472 ComplexNDArray::cummin (Array<octave_idx_type>& idx_arg, int dim) const
473 {
474 return do_mx_cumminmax_op<Complex> (*this, idx_arg, dim, mx_inline_cummin);
475 }
476
477 NDArray
abs(void) const478 ComplexNDArray::abs (void) const
479 {
480 return do_mx_unary_map<double, Complex, std::abs> (*this);
481 }
482
483 boolNDArray
isnan(void) const484 ComplexNDArray::isnan (void) const
485 {
486 return do_mx_unary_map<bool, Complex, octave::math::isnan> (*this);
487 }
488
489 boolNDArray
isinf(void) const490 ComplexNDArray::isinf (void) const
491 {
492 return do_mx_unary_map<bool, Complex, octave::math::isinf> (*this);
493 }
494
495 boolNDArray
isfinite(void) const496 ComplexNDArray::isfinite (void) const
497 {
498 return do_mx_unary_map<bool, Complex, octave::math::isfinite> (*this);
499 }
500
501 ComplexNDArray
conj(const ComplexNDArray & a)502 conj (const ComplexNDArray& a)
503 {
504 return do_mx_unary_map<Complex, Complex, std::conj<double>> (a);
505 }
506
507 ComplexNDArray&
insert(const NDArray & a,octave_idx_type r,octave_idx_type c)508 ComplexNDArray::insert (const NDArray& a, octave_idx_type r, octave_idx_type c)
509 {
510 dim_vector a_dv = a.dims ();
511
512 int n = a_dv.ndims ();
513
514 if (n != dimensions.ndims ())
515 (*current_liboctave_error_handler)
516 ("Array<T>::insert: invalid indexing operation");
517
518 Array<octave_idx_type> a_ra_idx (dim_vector (a_dv.ndims (), 1), 0);
519
520 a_ra_idx.elem (0) = r;
521 a_ra_idx.elem (1) = c;
522
523 for (int i = 0; i < n; i++)
524 {
525 if (a_ra_idx(i) < 0 || (a_ra_idx(i) + a_dv(i)) > dimensions(i))
526 (*current_liboctave_error_handler)
527 ("Array<T>::insert: range error for insert");
528 }
529
530 a_ra_idx.elem (0) = 0;
531 a_ra_idx.elem (1) = 0;
532
533 octave_idx_type n_elt = a.numel ();
534
535 // IS make_unique () NECESSARY HERE?
536
537 for (octave_idx_type i = 0; i < n_elt; i++)
538 {
539 Array<octave_idx_type> ra_idx = a_ra_idx;
540
541 ra_idx.elem (0) = a_ra_idx(0) + r;
542 ra_idx.elem (1) = a_ra_idx(1) + c;
543
544 elem (ra_idx) = a.elem (a_ra_idx);
545
546 increment_index (a_ra_idx, a_dv);
547 }
548
549 return *this;
550 }
551
552 ComplexNDArray&
insert(const ComplexNDArray & a,octave_idx_type r,octave_idx_type c)553 ComplexNDArray::insert (const ComplexNDArray& a,
554 octave_idx_type r, octave_idx_type c)
555 {
556 Array<Complex>::insert (a, r, c);
557 return *this;
558 }
559
560 ComplexNDArray&
insert(const ComplexNDArray & a,const Array<octave_idx_type> & ra_idx)561 ComplexNDArray::insert (const ComplexNDArray& a,
562 const Array<octave_idx_type>& ra_idx)
563 {
564 Array<Complex>::insert (a, ra_idx);
565 return *this;
566 }
567
568 void
increment_index(Array<octave_idx_type> & ra_idx,const dim_vector & dimensions,int start_dimension)569 ComplexNDArray::increment_index (Array<octave_idx_type>& ra_idx,
570 const dim_vector& dimensions,
571 int start_dimension)
572 {
573 ::increment_index (ra_idx, dimensions, start_dimension);
574 }
575
576 octave_idx_type
compute_index(Array<octave_idx_type> & ra_idx,const dim_vector & dimensions)577 ComplexNDArray::compute_index (Array<octave_idx_type>& ra_idx,
578 const dim_vector& dimensions)
579 {
580 return ::compute_index (ra_idx, dimensions);
581 }
582
583 ComplexNDArray
diag(octave_idx_type k) const584 ComplexNDArray::diag (octave_idx_type k) const
585 {
586 return MArray<Complex>::diag (k);
587 }
588
589 ComplexNDArray
diag(octave_idx_type m,octave_idx_type n) const590 ComplexNDArray::diag (octave_idx_type m, octave_idx_type n) const
591 {
592 return MArray<Complex>::diag (m, n);
593 }
594
595 // This contains no information on the array structure !!!
596 std::ostream&
operator <<(std::ostream & os,const ComplexNDArray & a)597 operator << (std::ostream& os, const ComplexNDArray& a)
598 {
599 octave_idx_type nel = a.numel ();
600
601 for (octave_idx_type i = 0; i < nel; i++)
602 {
603 os << ' ';
604 octave_write_complex (os, a.elem (i));
605 os << "\n";
606 }
607 return os;
608 }
609
610 std::istream&
operator >>(std::istream & is,ComplexNDArray & a)611 operator >> (std::istream& is, ComplexNDArray& a)
612 {
613 octave_idx_type nel = a.numel ();
614
615 if (nel > 0)
616 {
617 Complex tmp;
618 for (octave_idx_type i = 0; i < nel; i++)
619 {
620 tmp = octave_read_value<Complex> (is);
621 if (is)
622 a.elem (i) = tmp;
623 else
624 return is;
625 }
626 }
627
628 return is;
629 }
630
MINMAX_FCNS(ComplexNDArray,Complex)631 MINMAX_FCNS (ComplexNDArray, Complex)
632
633 NDS_CMP_OPS (ComplexNDArray, Complex)
634 NDS_BOOL_OPS (ComplexNDArray, Complex)
635
636 SND_CMP_OPS (Complex, ComplexNDArray)
637 SND_BOOL_OPS (Complex, ComplexNDArray)
638
639 NDND_CMP_OPS (ComplexNDArray, ComplexNDArray)
640 NDND_BOOL_OPS (ComplexNDArray, ComplexNDArray)
641
642 ComplexNDArray& operator *= (ComplexNDArray& a, double s)
643 {
644 if (a.is_shared ())
645 a = a * s;
646 else
647 do_ms_inplace_op<Complex, double> (a, s, mx_inline_mul2);
648 return a;
649 }
650
operator /=(ComplexNDArray & a,double s)651 ComplexNDArray& operator /= (ComplexNDArray& a, double s)
652 {
653 if (a.is_shared ())
654 return a = a / s;
655 else
656 do_ms_inplace_op<Complex, double> (a, s, mx_inline_div2);
657 return a;
658 }
659
660 BSXFUN_STDOP_DEFS_MXLOOP (ComplexNDArray)
661 BSXFUN_STDREL_DEFS_MXLOOP (ComplexNDArray)
662
663 BSXFUN_OP_DEF_MXLOOP (pow, ComplexNDArray, mx_inline_pow)
664