00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #if defined (__GNUG__) && defined (USE_PRAGMA_INTERFACE_IMPLEMENTATION)
00025 #pragma implementation
00026 #endif
00027
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031
00032 #include <cfloat>
00033 #include <vector>
00034
00035 #include "Array-util.h"
00036 #include "CNDArray.h"
00037 #include "mx-base.h"
00038 #include "f77-fcn.h"
00039 #include "lo-ieee.h"
00040 #include "lo-mappers.h"
00041
00042 #if defined (HAVE_FFTW3)
00043 #include "oct-fftw.h"
00044 #else
00045 extern "C"
00046 {
00047
00048
00049
00050
00051
00052 F77_RET_T
00053 F77_FUNC (cffti, CFFTI) (const int&, Complex*);
00054
00055 F77_RET_T
00056 F77_FUNC (cfftf, CFFTF) (const int&, Complex*, Complex*);
00057
00058 F77_RET_T
00059 F77_FUNC (cfftb, CFFTB) (const int&, Complex*, Complex*);
00060 }
00061 #endif
00062
00063 #if defined (HAVE_FFTW3)
00064 ComplexNDArray
00065 ComplexNDArray::fourier (int dim) const
00066 {
00067 dim_vector dv = dims ();
00068
00069 if (dim > dv.length () || dim < 0)
00070 return ComplexNDArray ();
00071
00072 int stride = 1;
00073 int n = dv(dim);
00074
00075 for (int i = 0; i < dim; i++)
00076 stride *= dv(i);
00077
00078 int howmany = numel () / dv (dim);
00079 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
00080 int nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride);
00081 int dist = (stride == 1 ? n : 1);
00082
00083 const Complex *in (fortran_vec ());
00084 ComplexNDArray retval (dv);
00085 Complex *out (retval.fortran_vec ());
00086
00087
00088 for (int k = 0; k < nloop; k++)
00089 octave_fftw::fft (in + k * stride * n, out + k * stride * n,
00090 n, howmany, stride, dist);
00091
00092 return retval;
00093 }
00094
00095 ComplexNDArray
00096 ComplexNDArray::ifourier (int dim) const
00097 {
00098 dim_vector dv = dims ();
00099
00100 if (dim > dv.length () || dim < 0)
00101 return ComplexNDArray ();
00102
00103 int stride = 1;
00104 int n = dv(dim);
00105
00106 for (int i = 0; i < dim; i++)
00107 stride *= dv(i);
00108
00109 int howmany = numel () / dv (dim);
00110 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
00111 int nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride);
00112 int dist = (stride == 1 ? n : 1);
00113
00114 const Complex *in (fortran_vec ());
00115 ComplexNDArray retval (dv);
00116 Complex *out (retval.fortran_vec ());
00117
00118
00119 for (int k = 0; k < nloop; k++)
00120 octave_fftw::ifft (in + k * stride * n, out + k * stride * n,
00121 n, howmany, stride, dist);
00122
00123 return retval;
00124 }
00125
00126 ComplexNDArray
00127 ComplexNDArray::fourier2d (void) const
00128 {
00129 dim_vector dv = dims();
00130 if (dv.length () < 2)
00131 return ComplexNDArray ();
00132
00133 dim_vector dv2(dv(0), dv(1));
00134 const Complex *in = fortran_vec ();
00135 ComplexNDArray retval (dv);
00136 Complex *out = retval.fortran_vec ();
00137 int howmany = numel() / dv(0) / dv(1);
00138 int dist = dv(0) * dv(1);
00139
00140 for (int i=0; i < howmany; i++)
00141 octave_fftw::fftNd (in + i*dist, out + i*dist, 2, dv2);
00142
00143 return retval;
00144 }
00145
00146 ComplexNDArray
00147 ComplexNDArray::ifourier2d (void) const
00148 {
00149 dim_vector dv = dims();
00150 if (dv.length () < 2)
00151 return ComplexNDArray ();
00152
00153 dim_vector dv2(dv(0), dv(1));
00154 const Complex *in = fortran_vec ();
00155 ComplexNDArray retval (dv);
00156 Complex *out = retval.fortran_vec ();
00157 int howmany = numel() / dv(0) / dv(1);
00158 int dist = dv(0) * dv(1);
00159
00160 for (int i=0; i < howmany; i++)
00161 octave_fftw::ifftNd (in + i*dist, out + i*dist, 2, dv2);
00162
00163 return retval;
00164 }
00165
00166 ComplexNDArray
00167 ComplexNDArray::fourierNd (void) const
00168 {
00169 dim_vector dv = dims ();
00170 int rank = dv.length ();
00171
00172 const Complex *in (fortran_vec ());
00173 ComplexNDArray retval (dv);
00174 Complex *out (retval.fortran_vec ());
00175
00176 octave_fftw::fftNd (in, out, rank, dv);
00177
00178 return retval;
00179 }
00180
00181 ComplexNDArray
00182 ComplexNDArray::ifourierNd (void) const
00183 {
00184 dim_vector dv = dims ();
00185 int rank = dv.length ();
00186
00187 const Complex *in (fortran_vec ());
00188 ComplexNDArray retval (dv);
00189 Complex *out (retval.fortran_vec ());
00190
00191 octave_fftw::ifftNd (in, out, rank, dv);
00192
00193 return retval;
00194 }
00195
00196 #else
00197 ComplexNDArray
00198 ComplexNDArray::fourier (int dim) const
00199 {
00200 dim_vector dv = dims ();
00201
00202 if (dim > dv.length () || dim < 0)
00203 return ComplexNDArray ();
00204
00205 ComplexNDArray retval (dv);
00206 int npts = dv(dim);
00207 int nn = 4*npts+15;
00208 Array<Complex> wsave (nn);
00209 Complex *pwsave = wsave.fortran_vec ();
00210
00211 OCTAVE_LOCAL_BUFFER (Complex, tmp, npts);
00212
00213 int stride = 1;
00214
00215 for (int i = 0; i < dim; i++)
00216 stride *= dv(i);
00217
00218 int howmany = numel () / npts;
00219 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
00220 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
00221 int dist = (stride == 1 ? npts : 1);
00222
00223 F77_FUNC (cffti, CFFTI) (npts, pwsave);
00224
00225 for (int k = 0; k < nloop; k++)
00226 {
00227 for (int j = 0; j < howmany; j++)
00228 {
00229 OCTAVE_QUIT;
00230
00231 for (int i = 0; i < npts; i++)
00232 tmp[i] = elem((i + k*npts)*stride + j*dist);
00233
00234 F77_FUNC (cfftf, CFFTF) (npts, tmp, pwsave);
00235
00236 for (int i = 0; i < npts; i++)
00237 retval ((i + k*npts)*stride + j*dist) = tmp[i];
00238 }
00239 }
00240
00241 return retval;
00242 }
00243
00244 ComplexNDArray
00245 ComplexNDArray::ifourier (int dim) const
00246 {
00247 dim_vector dv = dims ();
00248
00249 if (dim > dv.length () || dim < 0)
00250 return ComplexNDArray ();
00251
00252 ComplexNDArray retval (dv);
00253 int npts = dv(dim);
00254 int nn = 4*npts+15;
00255 Array<Complex> wsave (nn);
00256 Complex *pwsave = wsave.fortran_vec ();
00257
00258 OCTAVE_LOCAL_BUFFER (Complex, tmp, npts);
00259
00260 int stride = 1;
00261
00262 for (int i = 0; i < dim; i++)
00263 stride *= dv(i);
00264
00265 int howmany = numel () / npts;
00266 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
00267 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
00268 int dist = (stride == 1 ? npts : 1);
00269
00270 F77_FUNC (cffti, CFFTI) (npts, pwsave);
00271
00272 for (int k = 0; k < nloop; k++)
00273 {
00274 for (int j = 0; j < howmany; j++)
00275 {
00276 OCTAVE_QUIT;
00277
00278 for (int i = 0; i < npts; i++)
00279 tmp[i] = elem((i + k*npts)*stride + j*dist);
00280
00281 F77_FUNC (cfftb, CFFTB) (npts, tmp, pwsave);
00282
00283 for (int i = 0; i < npts; i++)
00284 retval ((i + k*npts)*stride + j*dist) = tmp[i] /
00285 static_cast<double> (npts);
00286 }
00287 }
00288
00289 return retval;
00290 }
00291
00292 ComplexNDArray
00293 ComplexNDArray::fourier2d (void) const
00294 {
00295 dim_vector dv = dims ();
00296 dim_vector dv2 (dv(0), dv(1));
00297 int rank = 2;
00298 ComplexNDArray retval (*this);
00299 int stride = 1;
00300
00301 for (int i = 0; i < rank; i++)
00302 {
00303 int npts = dv2(i);
00304 int nn = 4*npts+15;
00305 Array<Complex> wsave (nn);
00306 Complex *pwsave = wsave.fortran_vec ();
00307 Array<Complex> row (npts);
00308 Complex *prow = row.fortran_vec ();
00309
00310 int howmany = numel () / npts;
00311 howmany = (stride == 1 ? howmany :
00312 (howmany > stride ? stride : howmany));
00313 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
00314 int dist = (stride == 1 ? npts : 1);
00315
00316 F77_FUNC (cffti, CFFTI) (npts, pwsave);
00317
00318 for (int k = 0; k < nloop; k++)
00319 {
00320 for (int j = 0; j < howmany; j++)
00321 {
00322 OCTAVE_QUIT;
00323
00324 for (int l = 0; l < npts; l++)
00325 prow[l] = retval ((l + k*npts)*stride + j*dist);
00326
00327 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
00328
00329 for (int l = 0; l < npts; l++)
00330 retval ((l + k*npts)*stride + j*dist) = prow[l];
00331 }
00332 }
00333
00334 stride *= dv2(i);
00335 }
00336
00337 return retval;
00338 }
00339
00340 ComplexNDArray
00341 ComplexNDArray::ifourier2d (void) const
00342 {
00343 dim_vector dv = dims();
00344 dim_vector dv2 (dv(0), dv(1));
00345 int rank = 2;
00346 ComplexNDArray retval (*this);
00347 int stride = 1;
00348
00349 for (int i = 0; i < rank; i++)
00350 {
00351 int npts = dv2(i);
00352 int nn = 4*npts+15;
00353 Array<Complex> wsave (nn);
00354 Complex *pwsave = wsave.fortran_vec ();
00355 Array<Complex> row (npts);
00356 Complex *prow = row.fortran_vec ();
00357
00358 int howmany = numel () / npts;
00359 howmany = (stride == 1 ? howmany :
00360 (howmany > stride ? stride : howmany));
00361 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
00362 int dist = (stride == 1 ? npts : 1);
00363
00364 F77_FUNC (cffti, CFFTI) (npts, pwsave);
00365
00366 for (int k = 0; k < nloop; k++)
00367 {
00368 for (int j = 0; j < howmany; j++)
00369 {
00370 OCTAVE_QUIT;
00371
00372 for (int l = 0; l < npts; l++)
00373 prow[l] = retval ((l + k*npts)*stride + j*dist);
00374
00375 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
00376
00377 for (int l = 0; l < npts; l++)
00378 retval ((l + k*npts)*stride + j*dist) = prow[l] /
00379 static_cast<double> (npts);
00380 }
00381 }
00382
00383 stride *= dv2(i);
00384 }
00385
00386 return retval;
00387 }
00388
00389 ComplexNDArray
00390 ComplexNDArray::fourierNd (void) const
00391 {
00392 dim_vector dv = dims ();
00393 int rank = dv.length ();
00394 ComplexNDArray retval (*this);
00395 int stride = 1;
00396
00397 for (int i = 0; i < rank; i++)
00398 {
00399 int npts = dv(i);
00400 int nn = 4*npts+15;
00401 Array<Complex> wsave (nn);
00402 Complex *pwsave = wsave.fortran_vec ();
00403 Array<Complex> row (npts);
00404 Complex *prow = row.fortran_vec ();
00405
00406 int howmany = numel () / npts;
00407 howmany = (stride == 1 ? howmany :
00408 (howmany > stride ? stride : howmany));
00409 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
00410 int dist = (stride == 1 ? npts : 1);
00411
00412 F77_FUNC (cffti, CFFTI) (npts, pwsave);
00413
00414 for (int k = 0; k < nloop; k++)
00415 {
00416 for (int j = 0; j < howmany; j++)
00417 {
00418 OCTAVE_QUIT;
00419
00420 for (int l = 0; l < npts; l++)
00421 prow[l] = retval ((l + k*npts)*stride + j*dist);
00422
00423 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
00424
00425 for (int l = 0; l < npts; l++)
00426 retval ((l + k*npts)*stride + j*dist) = prow[l];
00427 }
00428 }
00429
00430 stride *= dv(i);
00431 }
00432
00433 return retval;
00434 }
00435
00436 ComplexNDArray
00437 ComplexNDArray::ifourierNd (void) const
00438 {
00439 dim_vector dv = dims ();
00440 int rank = dv.length ();
00441 ComplexNDArray retval (*this);
00442 int stride = 1;
00443
00444 for (int i = 0; i < rank; i++)
00445 {
00446 int npts = dv(i);
00447 int nn = 4*npts+15;
00448 Array<Complex> wsave (nn);
00449 Complex *pwsave = wsave.fortran_vec ();
00450 Array<Complex> row (npts);
00451 Complex *prow = row.fortran_vec ();
00452
00453 int howmany = numel () / npts;
00454 howmany = (stride == 1 ? howmany :
00455 (howmany > stride ? stride : howmany));
00456 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
00457 int dist = (stride == 1 ? npts : 1);
00458
00459 F77_FUNC (cffti, CFFTI) (npts, pwsave);
00460
00461 for (int k = 0; k < nloop; k++)
00462 {
00463 for (int j = 0; j < howmany; j++)
00464 {
00465 OCTAVE_QUIT;
00466
00467 for (int l = 0; l < npts; l++)
00468 prow[l] = retval ((l + k*npts)*stride + j*dist);
00469
00470 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
00471
00472 for (int l = 0; l < npts; l++)
00473 retval ((l + k*npts)*stride + j*dist) = prow[l] /
00474 static_cast<double> (npts);
00475 }
00476 }
00477
00478 stride *= dv(i);
00479 }
00480
00481 return retval;
00482 }
00483
00484 #endif
00485
00486
00487
00488 boolNDArray
00489 ComplexNDArray::operator ! (void) const
00490 {
00491 boolNDArray b (dims ());
00492
00493 for (int i = 0; i < length (); i++)
00494 b.elem (i) = elem (i) != 0.0;
00495
00496 return b;
00497 }
00498
00499
00500
00501 bool
00502 ComplexNDArray::any_element_is_inf_or_nan (void) const
00503 {
00504 int nel = nelem ();
00505
00506 for (int i = 0; i < nel; i++)
00507 {
00508 Complex val = elem (i);
00509 if (xisinf (val) || xisnan (val))
00510 return true;
00511 }
00512 return false;
00513 }
00514
00515
00516
00517 bool
00518 ComplexNDArray::all_elements_are_real (void) const
00519 {
00520 int nel = nelem ();
00521
00522 for (int i = 0; i < nel; i++)
00523 {
00524 double ip = imag (elem (i));
00525
00526 if (ip != 0.0 || lo_ieee_signbit (ip))
00527 return false;
00528 }
00529
00530 return true;
00531 }
00532
00533
00534
00535
00536
00537 bool
00538 ComplexNDArray::all_integers (double& max_val, double& min_val) const
00539 {
00540 int nel = nelem ();
00541
00542 if (nel > 0)
00543 {
00544 Complex val = elem (0);
00545
00546 double r_val = real (val);
00547 double i_val = imag (val);
00548
00549 max_val = r_val;
00550 min_val = r_val;
00551
00552 if (i_val > max_val)
00553 max_val = i_val;
00554
00555 if (i_val < max_val)
00556 min_val = i_val;
00557 }
00558 else
00559 return false;
00560
00561 for (int i = 0; i < nel; i++)
00562 {
00563 Complex val = elem (i);
00564
00565 double r_val = real (val);
00566 double i_val = imag (val);
00567
00568 if (r_val > max_val)
00569 max_val = r_val;
00570
00571 if (i_val > max_val)
00572 max_val = i_val;
00573
00574 if (r_val < min_val)
00575 min_val = r_val;
00576
00577 if (i_val < min_val)
00578 min_val = i_val;
00579
00580 if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val)
00581 return false;
00582 }
00583
00584 return true;
00585 }
00586
00587 bool
00588 ComplexNDArray::too_large_for_float (void) const
00589 {
00590 int nel = nelem ();
00591
00592 for (int i = 0; i < nel; i++)
00593 {
00594 Complex val = elem (i);
00595
00596 double r_val = real (val);
00597 double i_val = imag (val);
00598
00599 if (r_val > FLT_MAX
00600 || i_val > FLT_MAX
00601 || r_val < FLT_MIN
00602 || i_val < FLT_MIN)
00603 return true;
00604 }
00605
00606 return false;
00607 }
00608
00609 boolNDArray
00610 ComplexNDArray::all (int dim) const
00611 {
00612 MX_ND_ANY_ALL_REDUCTION
00613 (MX_ND_ALL_EVAL (elem (iter_idx) == Complex (0, 0)), true);
00614 }
00615
00616 boolNDArray
00617 ComplexNDArray::any (int dim) const
00618 {
00619 MX_ND_ANY_ALL_REDUCTION
00620 (MX_ND_ANY_EVAL (elem (iter_idx) != Complex (0, 0)), false);
00621 }
00622
00623 ComplexNDArray
00624 ComplexNDArray::cumprod (int dim) const
00625 {
00626 MX_ND_CUMULATIVE_OP (ComplexNDArray, Complex, Complex (1, 0), *);
00627 }
00628
00629 ComplexNDArray
00630 ComplexNDArray::cumsum (int dim) const
00631 {
00632 MX_ND_CUMULATIVE_OP (ComplexNDArray, Complex, Complex (0, 0), +);
00633 }
00634
00635 ComplexNDArray
00636 ComplexNDArray::prod (int dim) const
00637 {
00638 MX_ND_COMPLEX_OP_REDUCTION (*= elem (iter_idx), Complex (1, 0));
00639 }
00640
00641 ComplexNDArray
00642 ComplexNDArray::sumsq (int dim) const
00643 {
00644 MX_ND_COMPLEX_OP_REDUCTION
00645 (+= imag (elem (iter_idx))
00646 ? elem (iter_idx) * conj (elem (iter_idx))
00647 : std::pow (elem (iter_idx), 2), Complex (0, 0));
00648 }
00649
00650 ComplexNDArray
00651 ComplexNDArray::sum (int dim) const
00652 {
00653 MX_ND_COMPLEX_OP_REDUCTION (+= elem (iter_idx), Complex (0, 0));
00654 }
00655
00656 ComplexNDArray
00657 ComplexNDArray::concat (const ComplexNDArray& rb, const Array<int>& ra_idx)
00658 {
00659 if (rb.numel () > 0)
00660 insert (rb, ra_idx);
00661 return *this;
00662 }
00663
00664 ComplexNDArray
00665 ComplexNDArray::concat (const NDArray& rb, const Array<int>& ra_idx)
00666 {
00667 ComplexNDArray tmp (rb);
00668 if (rb.numel () > 0)
00669 insert (tmp, ra_idx);
00670 return *this;
00671 }
00672
00673 ComplexNDArray
00674 concat (NDArray& ra, ComplexNDArray& rb, const Array<int>& ra_idx)
00675 {
00676 ComplexNDArray retval (ra);
00677 if (rb.numel () > 0)
00678 retval.insert (rb, ra_idx);
00679 return retval;
00680 }
00681
00682 static const Complex Complex_NaN_result (octave_NaN, octave_NaN);
00683
00684 ComplexNDArray
00685 ComplexNDArray::max (int dim) const
00686 {
00687 ArrayN<int> dummy_idx;
00688 return max (dummy_idx, dim);
00689 }
00690
00691 ComplexNDArray
00692 ComplexNDArray::max (ArrayN<int>& idx_arg, int dim) const
00693 {
00694 dim_vector dv = dims ();
00695 dim_vector dr = dims ();
00696
00697 if (dv.numel () == 0 || dim > dv.length () || dim < 0)
00698 return ComplexNDArray ();
00699
00700 dr(dim) = 1;
00701
00702 ComplexNDArray result (dr);
00703 idx_arg.resize (dr);
00704
00705 int x_stride = 1;
00706 int x_len = dv(dim);
00707 for (int i = 0; i < dim; i++)
00708 x_stride *= dv(i);
00709
00710 for (int i = 0; i < dr.numel (); i++)
00711 {
00712 int x_offset;
00713 if (x_stride == 1)
00714 x_offset = i * x_len;
00715 else
00716 {
00717 int x_offset2 = 0;
00718 x_offset = i;
00719 while (x_offset >= x_stride)
00720 {
00721 x_offset -= x_stride;
00722 x_offset2++;
00723 }
00724 x_offset += x_offset2 * x_stride * x_len;
00725 }
00726
00727 int idx_j;
00728
00729 Complex tmp_max;
00730
00731 double abs_max = octave_NaN;
00732
00733 for (idx_j = 0; idx_j < x_len; idx_j++)
00734 {
00735 tmp_max = elem (idx_j * x_stride + x_offset);
00736
00737 if (! octave_is_NaN_or_NA (tmp_max))
00738 {
00739 abs_max = ::abs(tmp_max);
00740 break;
00741 }
00742 }
00743
00744 for (int j = idx_j+1; j < x_len; j++)
00745 {
00746 Complex tmp = elem (j * x_stride + x_offset);
00747
00748 if (octave_is_NaN_or_NA (tmp))
00749 continue;
00750
00751 double abs_tmp = ::abs (tmp);
00752
00753 if (abs_tmp > abs_max)
00754 {
00755 idx_j = j;
00756 tmp_max = tmp;
00757 abs_max = abs_tmp;
00758 }
00759 }
00760
00761 if (octave_is_NaN_or_NA (tmp_max))
00762 {
00763 result.elem (i) = Complex_NaN_result;
00764 idx_arg.elem (i) = 0;
00765 }
00766 else
00767 {
00768 result.elem (i) = tmp_max;
00769 idx_arg.elem (i) = idx_j;
00770 }
00771 }
00772
00773 return result;
00774 }
00775
00776 ComplexNDArray
00777 ComplexNDArray::min (int dim) const
00778 {
00779 ArrayN<int> dummy_idx;
00780 return min (dummy_idx, dim);
00781 }
00782
00783 ComplexNDArray
00784 ComplexNDArray::min (ArrayN<int>& idx_arg, int dim) const
00785 {
00786 dim_vector dv = dims ();
00787 dim_vector dr = dims ();
00788
00789 if (dv.numel () == 0 || dim > dv.length () || dim < 0)
00790 return ComplexNDArray ();
00791
00792 dr(dim) = 1;
00793
00794 ComplexNDArray result (dr);
00795 idx_arg.resize (dr);
00796
00797 int x_stride = 1;
00798 int x_len = dv(dim);
00799 for (int i = 0; i < dim; i++)
00800 x_stride *= dv(i);
00801
00802 for (int i = 0; i < dr.numel (); i++)
00803 {
00804 int x_offset;
00805 if (x_stride == 1)
00806 x_offset = i * x_len;
00807 else
00808 {
00809 int x_offset2 = 0;
00810 x_offset = i;
00811 while (x_offset >= x_stride)
00812 {
00813 x_offset -= x_stride;
00814 x_offset2++;
00815 }
00816 x_offset += x_offset2 * x_stride * x_len;
00817 }
00818
00819 int idx_j;
00820
00821 Complex tmp_min;
00822
00823 double abs_min = octave_NaN;
00824
00825 for (idx_j = 0; idx_j < x_len; idx_j++)
00826 {
00827 tmp_min = elem (idx_j * x_stride + x_offset);
00828
00829 if (! octave_is_NaN_or_NA (tmp_min))
00830 {
00831 abs_min = ::abs(tmp_min);
00832 break;
00833 }
00834 }
00835
00836 for (int j = idx_j+1; j < x_len; j++)
00837 {
00838 Complex tmp = elem (j * x_stride + x_offset);
00839
00840 if (octave_is_NaN_or_NA (tmp))
00841 continue;
00842
00843 double abs_tmp = ::abs (tmp);
00844
00845 if (abs_tmp < abs_min)
00846 {
00847 idx_j = j;
00848 tmp_min = tmp;
00849 abs_min = abs_tmp;
00850 }
00851 }
00852
00853 if (octave_is_NaN_or_NA (tmp_min))
00854 {
00855 result.elem (i) = Complex_NaN_result;
00856 idx_arg.elem (i) = 0;
00857 }
00858 else
00859 {
00860 result.elem (i) = tmp_min;
00861 idx_arg.elem (i) = idx_j;
00862 }
00863 }
00864
00865 return result;
00866 }
00867
00868 NDArray
00869 ComplexNDArray::abs (void) const
00870 {
00871 NDArray retval (dims ());
00872
00873 int nel = nelem ();
00874
00875 for (int i = 0; i < nel; i++)
00876 retval(i) = ::abs (elem (i));
00877
00878 return retval;
00879 }
00880
00881 ComplexNDArray&
00882 ComplexNDArray::insert (const NDArray& a, int r, int c)
00883 {
00884 dim_vector a_dv = a.dims ();
00885
00886 int n = a_dv.length ();
00887
00888 if (n == dimensions.length ())
00889 {
00890 Array<int> a_ra_idx (a_dv.length (), 0);
00891
00892 a_ra_idx.elem (0) = r;
00893 a_ra_idx.elem (1) = c;
00894
00895 for (int i = 0; i < n; i++)
00896 {
00897 if (a_ra_idx (i) < 0 || (a_ra_idx (i) + a_dv (i)) > dimensions (i))
00898 {
00899 (*current_liboctave_error_handler)
00900 ("Array<T>::insert: range error for insert");
00901 return *this;
00902 }
00903 }
00904
00905 a_ra_idx.elem (0) = 0;
00906 a_ra_idx.elem (1) = 0;
00907
00908 int n_elt = a.numel ();
00909
00910
00911
00912 for (int i = 0; i < n_elt; i++)
00913 {
00914 Array<int> ra_idx = a_ra_idx;
00915
00916 ra_idx.elem (0) = a_ra_idx (0) + r;
00917 ra_idx.elem (1) = a_ra_idx (1) + c;
00918
00919 elem (ra_idx) = a.elem (a_ra_idx);
00920
00921 increment_index (a_ra_idx, a_dv);
00922 }
00923 }
00924 else
00925 (*current_liboctave_error_handler)
00926 ("Array<T>::insert: invalid indexing operation");
00927
00928 return *this;
00929 }
00930
00931 ComplexNDArray&
00932 ComplexNDArray::insert (const ComplexNDArray& a, int r, int c)
00933 {
00934 Array<Complex>::insert (a, r, c);
00935 return *this;
00936 }
00937
00938 ComplexNDArray&
00939 ComplexNDArray::insert (const ComplexNDArray& a, const Array<int>& ra_idx)
00940 {
00941 Array<Complex>::insert (a, ra_idx);
00942 return *this;
00943 }
00944
00945 ComplexMatrix
00946 ComplexNDArray::matrix_value (void) const
00947 {
00948 ComplexMatrix retval;
00949
00950 int nd = ndims ();
00951
00952 switch (nd)
00953 {
00954 case 1:
00955 retval = ComplexMatrix (Array2<Complex> (*this, dimensions(0), 1));
00956 break;
00957
00958 case 2:
00959 retval = ComplexMatrix (Array2<Complex> (*this, dimensions(0),
00960 dimensions(1)));
00961 break;
00962
00963 default:
00964 (*current_liboctave_error_handler)
00965 ("invalid conversion of ComplexNDArray to ComplexMatrix");
00966 break;
00967 }
00968
00969 return retval;
00970 }
00971
00972 void
00973 ComplexNDArray::increment_index (Array<int>& ra_idx,
00974 const dim_vector& dimensions,
00975 int start_dimension)
00976 {
00977 ::increment_index (ra_idx, dimensions, start_dimension);
00978 }
00979
00980 int
00981 ComplexNDArray::compute_index (Array<int>& ra_idx,
00982 const dim_vector& dimensions)
00983 {
00984 return ::compute_index (ra_idx, dimensions);
00985 }
00986
00987
00988
00989 std::ostream&
00990 operator << (std::ostream& os, const ComplexNDArray& a)
00991 {
00992 int nel = a.nelem ();
00993
00994 for (int i = 0; i < nel; i++)
00995 {
00996 os << " ";
00997 octave_write_complex (os, a.elem (i));
00998 os << "\n";
00999 }
01000 return os;
01001 }
01002
01003 std::istream&
01004 operator >> (std::istream& is, ComplexNDArray& a)
01005 {
01006 int nel = a.nelem ();
01007
01008 if (nel < 1 )
01009 is.clear (std::ios::badbit);
01010 else
01011 {
01012 Complex tmp;
01013 for (int i = 0; i < nel; i++)
01014 {
01015 tmp = octave_read_complex (is);
01016 if (is)
01017 a.elem (i) = tmp;
01018 else
01019 goto done;
01020 }
01021 }
01022
01023 done:
01024
01025 return is;
01026 }
01027
01028
01029
01030
01031 #define EMPTY_RETURN_CHECK(T) \
01032 if (nel == 0) \
01033 return T (dv);
01034
01035 ComplexNDArray
01036 min (const Complex& c, const ComplexNDArray& m)
01037 {
01038 dim_vector dv = m.dims ();
01039 int nel = dv.numel ();
01040
01041 EMPTY_RETURN_CHECK (ComplexNDArray);
01042
01043 ComplexNDArray result (dv);
01044
01045 for (int i = 0; i < nel; i++)
01046 {
01047 OCTAVE_QUIT;
01048 result (i) = xmin (c, m (i));
01049 }
01050
01051 return result;
01052 }
01053
01054 ComplexNDArray
01055 min (const ComplexNDArray& m, const Complex& c)
01056 {
01057 dim_vector dv = m.dims ();
01058 int nel = dv.numel ();
01059
01060 EMPTY_RETURN_CHECK (ComplexNDArray);
01061
01062 ComplexNDArray result (dv);
01063
01064 for (int i = 0; i < nel; i++)
01065 {
01066 OCTAVE_QUIT;
01067 result (i) = xmin (c, m (i));
01068 }
01069
01070 return result;
01071 }
01072
01073 ComplexNDArray
01074 min (const ComplexNDArray& a, const ComplexNDArray& b)
01075 {
01076 dim_vector dv = a.dims ();
01077 int nel = dv.numel ();
01078
01079 if (dv != b.dims ())
01080 {
01081 (*current_liboctave_error_handler)
01082 ("two-arg min expecting args of same size");
01083 return ComplexNDArray ();
01084 }
01085
01086 EMPTY_RETURN_CHECK (ComplexNDArray);
01087
01088 ComplexNDArray result (dv);
01089
01090 for (int i = 0; i < nel; i++)
01091 {
01092 OCTAVE_QUIT;
01093 result (i) = xmin (a (i), b (i));
01094 }
01095
01096 return result;
01097 }
01098
01099 ComplexNDArray
01100 max (const Complex& c, const ComplexNDArray& m)
01101 {
01102 dim_vector dv = m.dims ();
01103 int nel = dv.numel ();
01104
01105 EMPTY_RETURN_CHECK (ComplexNDArray);
01106
01107 ComplexNDArray result (dv);
01108
01109 for (int i = 0; i < nel; i++)
01110 {
01111 OCTAVE_QUIT;
01112 result (i) = xmax (c, m (i));
01113 }
01114
01115 return result;
01116 }
01117
01118 ComplexNDArray
01119 max (const ComplexNDArray& m, const Complex& c)
01120 {
01121 dim_vector dv = m.dims ();
01122 int nel = dv.numel ();
01123
01124 EMPTY_RETURN_CHECK (ComplexNDArray);
01125
01126 ComplexNDArray result (dv);
01127
01128 for (int i = 0; i < nel; i++)
01129 {
01130 OCTAVE_QUIT;
01131 result (i) = xmax (c, m (i));
01132 }
01133
01134 return result;
01135 }
01136
01137 ComplexNDArray
01138 max (const ComplexNDArray& a, const ComplexNDArray& b)
01139 {
01140 dim_vector dv = a.dims ();
01141 int nel = dv.numel ();
01142
01143 if (dv != b.dims ())
01144 {
01145 (*current_liboctave_error_handler)
01146 ("two-arg max expecting args of same size");
01147 return ComplexNDArray ();
01148 }
01149
01150 EMPTY_RETURN_CHECK (ComplexNDArray);
01151
01152 ComplexNDArray result (dv);
01153
01154 for (int i = 0; i < nel; i++)
01155 {
01156 OCTAVE_QUIT;
01157 result (i) = xmax (a (i), b (i));
01158 }
01159
01160 return result;
01161 }
01162
01163 NDS_CMP_OPS(ComplexNDArray, real, Complex, real)
01164 NDS_BOOL_OPS(ComplexNDArray, Complex, 0.0)
01165
01166 SND_CMP_OPS(Complex, real, ComplexNDArray, real)
01167 SND_BOOL_OPS(Complex, ComplexNDArray, 0.0)
01168
01169 NDND_CMP_OPS(ComplexNDArray, real, ComplexNDArray, real)
01170 NDND_BOOL_OPS(ComplexNDArray, ComplexNDArray, 0.0)
01171
01172
01173
01174
01175
01176