Mercurial > octave-nkf
annotate liboctave/dNDArray.cc @ 8751:9f7ce4bf7650
optimize min/max functions
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Mon, 16 Feb 2009 08:52:00 +0100 |
parents | 1bd918cfb6e2 |
children | b756ce0002db |
rev | line source |
---|---|
4513 | 1 // N-D Array manipulations. |
4511 | 2 /* |
3 | |
7017 | 4 Copyright (C) 1996, 1997, 2003, 2004, 2005, 2006, 2007 John W. Eaton |
4511 | 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 | |
7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
11 option) any later version. | |
4511 | 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 | |
7016 | 19 along with Octave; see the file COPYING. If not, see |
20 <http://www.gnu.org/licenses/>. | |
4511 | 21 |
22 */ | |
23 | |
24 #ifdef HAVE_CONFIG_H | |
25 #include <config.h> | |
26 #endif | |
27 | |
4634 | 28 #include <cfloat> |
5164 | 29 |
4780 | 30 #include <vector> |
4634 | 31 |
4588 | 32 #include "Array-util.h" |
4513 | 33 #include "dNDArray.h" |
7503
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
34 #include "functor.h" |
4511 | 35 #include "mx-base.h" |
4773 | 36 #include "f77-fcn.h" |
4513 | 37 #include "lo-error.h" |
4511 | 38 #include "lo-ieee.h" |
4634 | 39 #include "lo-mappers.h" |
8377
25bc2d31e1bf
improve OCTAVE_LOCAL_BUFFER
Jaroslav Hajek <highegg@gmail.com>
parents:
7941
diff
changeset
|
40 #include "oct-locbuf.h" |
4511 | 41 |
4773 | 42 #if defined (HAVE_FFTW3) |
43 #include "oct-fftw.h" | |
7941
f8cab9eeb128
Fix NDArray compilation/export
John W. Eaton <jwe@octave.org>
parents:
7922
diff
changeset
|
44 #endif |
4773 | 45 |
7919
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
46 NDArray::NDArray (const Array<octave_idx_type>& a, bool zero_based, |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
47 bool negative_to_nan) |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
48 { |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
49 const octave_idx_type *pa = a.fortran_vec (); |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
50 resize (a.dims ()); |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
51 double *ptmp = fortran_vec (); |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
52 if (negative_to_nan) |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
53 { |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
54 double nan_val = lo_ieee_nan_value (); |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
55 |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
56 if (zero_based) |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
57 for (octave_idx_type i = 0; i < a.numel (); i++) |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
58 { |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
59 double val = static_cast<double> |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
60 (pa[i] + static_cast<octave_idx_type> (1)); |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
61 if (val <= 0) |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
62 ptmp[i] = nan_val; |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
63 else |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
64 ptmp[i] = val; |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
65 } |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
66 else |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
67 for (octave_idx_type i = 0; i < a.numel (); i++) |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
68 { |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
69 double val = static_cast<double> (pa[i]); |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
70 if (val <= 0) |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
71 ptmp[i] = nan_val; |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
72 else |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
73 ptmp[i] = val; |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
74 } |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
75 } |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
76 else |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
77 { |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
78 if (zero_based) |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
79 for (octave_idx_type i = 0; i < a.numel (); i++) |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
80 ptmp[i] = static_cast<double> |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
81 (pa[i] + static_cast<octave_idx_type> (1)); |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
82 else |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
83 for (octave_idx_type i = 0; i < a.numel (); i++) |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
84 ptmp[i] = static_cast<double> (pa[i]); |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
85 } |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
86 } |
9d080df0c843
new NDArray constructor for ArrayN<octave_idx_type>
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
87 |
7941
f8cab9eeb128
Fix NDArray compilation/export
John W. Eaton <jwe@octave.org>
parents:
7922
diff
changeset
|
88 #if defined (HAVE_FFTW3) |
f8cab9eeb128
Fix NDArray compilation/export
John W. Eaton <jwe@octave.org>
parents:
7922
diff
changeset
|
89 |
4773 | 90 ComplexNDArray |
91 NDArray::fourier (int dim) const | |
92 { | |
93 dim_vector dv = dims (); | |
94 | |
95 if (dim > dv.length () || dim < 0) | |
96 return ComplexNDArray (); | |
97 | |
5275 | 98 octave_idx_type stride = 1; |
99 octave_idx_type n = dv(dim); | |
4773 | 100 |
101 for (int i = 0; i < dim; i++) | |
102 stride *= dv(i); | |
103 | |
5275 | 104 octave_idx_type howmany = numel () / dv (dim); |
4773 | 105 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany)); |
5275 | 106 octave_idx_type nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride); |
107 octave_idx_type dist = (stride == 1 ? n : 1); | |
4773 | 108 |
109 const double *in (fortran_vec ()); | |
110 ComplexNDArray retval (dv); | |
111 Complex *out (retval.fortran_vec ()); | |
112 | |
113 // Need to be careful here about the distance between fft's | |
5275 | 114 for (octave_idx_type k = 0; k < nloop; k++) |
4773 | 115 octave_fftw::fft (in + k * stride * n, out + k * stride * n, |
116 n, howmany, stride, dist); | |
117 | |
118 return retval; | |
119 } | |
120 | |
121 ComplexNDArray | |
4816 | 122 NDArray::ifourier (int dim) const |
4773 | 123 { |
124 dim_vector dv = dims (); | |
125 | |
126 if (dim > dv.length () || dim < 0) | |
127 return ComplexNDArray (); | |
128 | |
5275 | 129 octave_idx_type stride = 1; |
130 octave_idx_type n = dv(dim); | |
4773 | 131 |
132 for (int i = 0; i < dim; i++) | |
133 stride *= dv(i); | |
134 | |
5275 | 135 octave_idx_type howmany = numel () / dv (dim); |
4773 | 136 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany)); |
5275 | 137 octave_idx_type nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride); |
138 octave_idx_type dist = (stride == 1 ? n : 1); | |
4773 | 139 |
140 ComplexNDArray retval (*this); | |
141 Complex *out (retval.fortran_vec ()); | |
142 | |
143 // Need to be careful here about the distance between fft's | |
5275 | 144 for (octave_idx_type k = 0; k < nloop; k++) |
4773 | 145 octave_fftw::ifft (out + k * stride * n, out + k * stride * n, |
146 n, howmany, stride, dist); | |
147 | |
148 return retval; | |
149 } | |
150 | |
151 ComplexNDArray | |
152 NDArray::fourier2d (void) const | |
153 { | |
154 dim_vector dv = dims(); | |
155 if (dv.length () < 2) | |
156 return ComplexNDArray (); | |
157 | |
158 dim_vector dv2(dv(0), dv(1)); | |
159 const double *in = fortran_vec (); | |
160 ComplexNDArray retval (dv); | |
161 Complex *out = retval.fortran_vec (); | |
5275 | 162 octave_idx_type howmany = numel() / dv(0) / dv(1); |
163 octave_idx_type dist = dv(0) * dv(1); | |
4773 | 164 |
5275 | 165 for (octave_idx_type i=0; i < howmany; i++) |
4773 | 166 octave_fftw::fftNd (in + i*dist, out + i*dist, 2, dv2); |
167 | |
168 return retval; | |
169 } | |
170 | |
171 ComplexNDArray | |
172 NDArray::ifourier2d (void) const | |
173 { | |
174 dim_vector dv = dims(); | |
175 if (dv.length () < 2) | |
176 return ComplexNDArray (); | |
177 | |
178 dim_vector dv2(dv(0), dv(1)); | |
179 ComplexNDArray retval (*this); | |
180 Complex *out = retval.fortran_vec (); | |
5275 | 181 octave_idx_type howmany = numel() / dv(0) / dv(1); |
182 octave_idx_type dist = dv(0) * dv(1); | |
4773 | 183 |
5275 | 184 for (octave_idx_type i=0; i < howmany; i++) |
4773 | 185 octave_fftw::ifftNd (out + i*dist, out + i*dist, 2, dv2); |
186 | |
187 return retval; | |
188 } | |
189 | |
190 ComplexNDArray | |
191 NDArray::fourierNd (void) const | |
192 { | |
193 dim_vector dv = dims (); | |
194 int rank = dv.length (); | |
195 | |
196 const double *in (fortran_vec ()); | |
197 ComplexNDArray retval (dv); | |
198 Complex *out (retval.fortran_vec ()); | |
199 | |
200 octave_fftw::fftNd (in, out, rank, dv); | |
201 | |
202 return retval; | |
203 } | |
204 | |
205 ComplexNDArray | |
206 NDArray::ifourierNd (void) const | |
207 { | |
208 dim_vector dv = dims (); | |
209 int rank = dv.length (); | |
210 | |
211 ComplexNDArray tmp (*this); | |
212 Complex *in (tmp.fortran_vec ()); | |
213 ComplexNDArray retval (dv); | |
214 Complex *out (retval.fortran_vec ()); | |
215 | |
216 octave_fftw::ifftNd (in, out, rank, dv); | |
217 | |
218 return retval; | |
219 } | |
220 | |
221 #else | |
222 | |
223 extern "C" | |
224 { | |
225 // Note that the original complex fft routines were not written for | |
226 // double complex arguments. They have been modified by adding an | |
227 // implicit double precision (a-h,o-z) statement at the beginning of | |
228 // each subroutine. | |
229 | |
230 F77_RET_T | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
231 F77_FUNC (zffti, ZFFTI) (const octave_idx_type&, Complex*); |
4773 | 232 |
233 F77_RET_T | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
234 F77_FUNC (zfftf, ZFFTF) (const octave_idx_type&, Complex*, Complex*); |
4773 | 235 |
236 F77_RET_T | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
237 F77_FUNC (zfftb, ZFFTB) (const octave_idx_type&, Complex*, Complex*); |
4773 | 238 } |
239 | |
240 ComplexNDArray | |
241 NDArray::fourier (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); | |
5275 | 249 octave_idx_type npts = dv(dim); |
250 octave_idx_type nn = 4*npts+15; | |
4773 | 251 Array<Complex> wsave (nn); |
252 Complex *pwsave = wsave.fortran_vec (); | |
253 | |
254 OCTAVE_LOCAL_BUFFER (Complex, tmp, npts); | |
255 | |
5275 | 256 octave_idx_type stride = 1; |
4773 | 257 |
258 for (int i = 0; i < dim; i++) | |
259 stride *= dv(i); | |
260 | |
5275 | 261 octave_idx_type howmany = numel () / npts; |
4773 | 262 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany)); |
5275 | 263 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride); |
264 octave_idx_type dist = (stride == 1 ? npts : 1); | |
4773 | 265 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
266 F77_FUNC (zffti, ZFFTI) (npts, pwsave); |
4773 | 267 |
5275 | 268 for (octave_idx_type k = 0; k < nloop; k++) |
4773 | 269 { |
5275 | 270 for (octave_idx_type j = 0; j < howmany; j++) |
4773 | 271 { |
272 OCTAVE_QUIT; | |
273 | |
5275 | 274 for (octave_idx_type i = 0; i < npts; i++) |
4773 | 275 tmp[i] = elem((i + k*npts)*stride + j*dist); |
276 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
277 F77_FUNC (zfftf, ZFFTF) (npts, tmp, pwsave); |
4773 | 278 |
5275 | 279 for (octave_idx_type i = 0; i < npts; i++) |
4773 | 280 retval ((i + k*npts)*stride + j*dist) = tmp[i]; |
281 } | |
282 } | |
283 | |
284 return retval; | |
285 } | |
286 | |
287 ComplexNDArray | |
288 NDArray::ifourier (int dim) const | |
289 { | |
290 dim_vector dv = dims (); | |
291 | |
292 if (dim > dv.length () || dim < 0) | |
293 return ComplexNDArray (); | |
294 | |
295 ComplexNDArray retval (dv); | |
5275 | 296 octave_idx_type npts = dv(dim); |
297 octave_idx_type nn = 4*npts+15; | |
4773 | 298 Array<Complex> wsave (nn); |
299 Complex *pwsave = wsave.fortran_vec (); | |
300 | |
301 OCTAVE_LOCAL_BUFFER (Complex, tmp, npts); | |
302 | |
5275 | 303 octave_idx_type stride = 1; |
4773 | 304 |
305 for (int i = 0; i < dim; i++) | |
306 stride *= dv(i); | |
307 | |
5275 | 308 octave_idx_type howmany = numel () / npts; |
4773 | 309 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany)); |
5275 | 310 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride); |
311 octave_idx_type dist = (stride == 1 ? npts : 1); | |
4773 | 312 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
313 F77_FUNC (zffti, ZFFTI) (npts, pwsave); |
4773 | 314 |
5275 | 315 for (octave_idx_type k = 0; k < nloop; k++) |
4773 | 316 { |
5275 | 317 for (octave_idx_type j = 0; j < howmany; j++) |
4773 | 318 { |
319 OCTAVE_QUIT; | |
320 | |
5275 | 321 for (octave_idx_type i = 0; i < npts; i++) |
4773 | 322 tmp[i] = elem((i + k*npts)*stride + j*dist); |
323 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
324 F77_FUNC (zfftb, ZFFTB) (npts, tmp, pwsave); |
4773 | 325 |
5275 | 326 for (octave_idx_type i = 0; i < npts; i++) |
4773 | 327 retval ((i + k*npts)*stride + j*dist) = tmp[i] / |
328 static_cast<double> (npts); | |
329 } | |
330 } | |
331 | |
332 return retval; | |
333 } | |
334 | |
335 ComplexNDArray | |
336 NDArray::fourier2d (void) const | |
337 { | |
338 dim_vector dv = dims(); | |
339 dim_vector dv2 (dv(0), dv(1)); | |
340 int rank = 2; | |
341 ComplexNDArray retval (*this); | |
5275 | 342 octave_idx_type stride = 1; |
4773 | 343 |
344 for (int i = 0; i < rank; i++) | |
345 { | |
5275 | 346 octave_idx_type npts = dv2(i); |
347 octave_idx_type nn = 4*npts+15; | |
4773 | 348 Array<Complex> wsave (nn); |
349 Complex *pwsave = wsave.fortran_vec (); | |
350 Array<Complex> row (npts); | |
351 Complex *prow = row.fortran_vec (); | |
352 | |
5275 | 353 octave_idx_type howmany = numel () / npts; |
4773 | 354 howmany = (stride == 1 ? howmany : |
355 (howmany > stride ? stride : howmany)); | |
5275 | 356 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride); |
357 octave_idx_type dist = (stride == 1 ? npts : 1); | |
4773 | 358 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
359 F77_FUNC (zffti, ZFFTI) (npts, pwsave); |
4773 | 360 |
5275 | 361 for (octave_idx_type k = 0; k < nloop; k++) |
4773 | 362 { |
5275 | 363 for (octave_idx_type j = 0; j < howmany; j++) |
4773 | 364 { |
365 OCTAVE_QUIT; | |
366 | |
5275 | 367 for (octave_idx_type l = 0; l < npts; l++) |
4773 | 368 prow[l] = retval ((l + k*npts)*stride + j*dist); |
369 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
370 F77_FUNC (zfftf, ZFFTF) (npts, prow, pwsave); |
4773 | 371 |
5275 | 372 for (octave_idx_type l = 0; l < npts; l++) |
4773 | 373 retval ((l + k*npts)*stride + j*dist) = prow[l]; |
374 } | |
375 } | |
376 | |
377 stride *= dv2(i); | |
378 } | |
379 | |
380 return retval; | |
381 } | |
382 | |
383 ComplexNDArray | |
384 NDArray::ifourier2d (void) const | |
385 { | |
386 dim_vector dv = dims(); | |
387 dim_vector dv2 (dv(0), dv(1)); | |
388 int rank = 2; | |
389 ComplexNDArray retval (*this); | |
5275 | 390 octave_idx_type stride = 1; |
4773 | 391 |
392 for (int i = 0; i < rank; i++) | |
393 { | |
5275 | 394 octave_idx_type npts = dv2(i); |
395 octave_idx_type nn = 4*npts+15; | |
4773 | 396 Array<Complex> wsave (nn); |
397 Complex *pwsave = wsave.fortran_vec (); | |
398 Array<Complex> row (npts); | |
399 Complex *prow = row.fortran_vec (); | |
400 | |
5275 | 401 octave_idx_type howmany = numel () / npts; |
4773 | 402 howmany = (stride == 1 ? howmany : |
403 (howmany > stride ? stride : howmany)); | |
5275 | 404 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride); |
405 octave_idx_type dist = (stride == 1 ? npts : 1); | |
4773 | 406 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
407 F77_FUNC (zffti, ZFFTI) (npts, pwsave); |
4773 | 408 |
5275 | 409 for (octave_idx_type k = 0; k < nloop; k++) |
4773 | 410 { |
5275 | 411 for (octave_idx_type j = 0; j < howmany; j++) |
4773 | 412 { |
413 OCTAVE_QUIT; | |
414 | |
5275 | 415 for (octave_idx_type l = 0; l < npts; l++) |
4773 | 416 prow[l] = retval ((l + k*npts)*stride + j*dist); |
417 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
418 F77_FUNC (zfftb, ZFFTB) (npts, prow, pwsave); |
4773 | 419 |
5275 | 420 for (octave_idx_type l = 0; l < npts; l++) |
4773 | 421 retval ((l + k*npts)*stride + j*dist) = prow[l] / |
422 static_cast<double> (npts); | |
423 } | |
424 } | |
425 | |
426 stride *= dv2(i); | |
427 } | |
428 | |
429 return retval; | |
430 } | |
431 | |
432 ComplexNDArray | |
433 NDArray::fourierNd (void) const | |
434 { | |
435 dim_vector dv = dims (); | |
436 int rank = dv.length (); | |
437 ComplexNDArray retval (*this); | |
5275 | 438 octave_idx_type stride = 1; |
4773 | 439 |
440 for (int i = 0; i < rank; i++) | |
441 { | |
5275 | 442 octave_idx_type npts = dv(i); |
443 octave_idx_type nn = 4*npts+15; | |
4773 | 444 Array<Complex> wsave (nn); |
445 Complex *pwsave = wsave.fortran_vec (); | |
446 Array<Complex> row (npts); | |
447 Complex *prow = row.fortran_vec (); | |
448 | |
5275 | 449 octave_idx_type howmany = numel () / npts; |
4773 | 450 howmany = (stride == 1 ? howmany : |
451 (howmany > stride ? stride : howmany)); | |
5275 | 452 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride); |
453 octave_idx_type dist = (stride == 1 ? npts : 1); | |
4773 | 454 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
455 F77_FUNC (zffti, ZFFTI) (npts, pwsave); |
4773 | 456 |
5275 | 457 for (octave_idx_type k = 0; k < nloop; k++) |
4773 | 458 { |
5275 | 459 for (octave_idx_type j = 0; j < howmany; j++) |
4773 | 460 { |
461 OCTAVE_QUIT; | |
462 | |
5275 | 463 for (octave_idx_type l = 0; l < npts; l++) |
4773 | 464 prow[l] = retval ((l + k*npts)*stride + j*dist); |
465 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
466 F77_FUNC (zfftf, ZFFTF) (npts, prow, pwsave); |
4773 | 467 |
5275 | 468 for (octave_idx_type l = 0; l < npts; l++) |
4773 | 469 retval ((l + k*npts)*stride + j*dist) = prow[l]; |
470 } | |
471 } | |
472 | |
473 stride *= dv(i); | |
474 } | |
475 | |
476 return retval; | |
477 } | |
478 | |
479 ComplexNDArray | |
480 NDArray::ifourierNd (void) const | |
481 { | |
482 dim_vector dv = dims (); | |
483 int rank = dv.length (); | |
484 ComplexNDArray retval (*this); | |
5275 | 485 octave_idx_type stride = 1; |
4773 | 486 |
487 for (int i = 0; i < rank; i++) | |
488 { | |
5275 | 489 octave_idx_type npts = dv(i); |
490 octave_idx_type nn = 4*npts+15; | |
4773 | 491 Array<Complex> wsave (nn); |
492 Complex *pwsave = wsave.fortran_vec (); | |
493 Array<Complex> row (npts); | |
494 Complex *prow = row.fortran_vec (); | |
495 | |
5275 | 496 octave_idx_type howmany = numel () / npts; |
4773 | 497 howmany = (stride == 1 ? howmany : |
498 (howmany > stride ? stride : howmany)); | |
5275 | 499 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride); |
500 octave_idx_type dist = (stride == 1 ? npts : 1); | |
4773 | 501 |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
502 F77_FUNC (zffti, ZFFTI) (npts, pwsave); |
4773 | 503 |
5275 | 504 for (octave_idx_type k = 0; k < nloop; k++) |
4773 | 505 { |
5275 | 506 for (octave_idx_type j = 0; j < howmany; j++) |
4773 | 507 { |
508 OCTAVE_QUIT; | |
509 | |
5275 | 510 for (octave_idx_type l = 0; l < npts; l++) |
4773 | 511 prow[l] = retval ((l + k*npts)*stride + j*dist); |
512 | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
513 F77_FUNC (zfftb, ZFFTB) (npts, prow, pwsave); |
4773 | 514 |
5275 | 515 for (octave_idx_type l = 0; l < npts; l++) |
4773 | 516 retval ((l + k*npts)*stride + j*dist) = prow[l] / |
517 static_cast<double> (npts); | |
518 } | |
519 } | |
520 | |
521 stride *= dv(i); | |
522 } | |
523 | |
524 return retval; | |
525 } | |
526 | |
527 #endif | |
528 | |
4543 | 529 // unary operations |
530 | |
531 boolNDArray | |
532 NDArray::operator ! (void) const | |
533 { | |
534 boolNDArray b (dims ()); | |
535 | |
5275 | 536 for (octave_idx_type i = 0; i < length (); i++) |
4543 | 537 b.elem (i) = ! elem (i); |
538 | |
539 return b; | |
540 } | |
541 | |
4634 | 542 bool |
543 NDArray::any_element_is_negative (bool neg_zero) const | |
544 { | |
5275 | 545 octave_idx_type nel = nelem (); |
4634 | 546 |
547 if (neg_zero) | |
548 { | |
5275 | 549 for (octave_idx_type i = 0; i < nel; i++) |
4634 | 550 if (lo_ieee_signbit (elem (i))) |
551 return true; | |
552 } | |
553 else | |
554 { | |
5275 | 555 for (octave_idx_type i = 0; i < nel; i++) |
4634 | 556 if (elem (i) < 0) |
557 return true; | |
558 } | |
559 | |
560 return false; | |
561 } | |
562 | |
7922
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7919
diff
changeset
|
563 bool |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7919
diff
changeset
|
564 NDArray::any_element_is_nan (void) const |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7919
diff
changeset
|
565 { |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7919
diff
changeset
|
566 octave_idx_type nel = nelem (); |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7919
diff
changeset
|
567 |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7919
diff
changeset
|
568 for (octave_idx_type i = 0; i < nel; i++) |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7919
diff
changeset
|
569 { |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7919
diff
changeset
|
570 double val = elem (i); |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7919
diff
changeset
|
571 if (xisnan (val)) |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7919
diff
changeset
|
572 return true; |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7919
diff
changeset
|
573 } |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7919
diff
changeset
|
574 |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7919
diff
changeset
|
575 return false; |
935be827eaf8
error for NaN values in & and | expressions
John W. Eaton <jwe@octave.org>
parents:
7919
diff
changeset
|
576 } |
4634 | 577 |
578 bool | |
579 NDArray::any_element_is_inf_or_nan (void) const | |
580 { | |
5275 | 581 octave_idx_type nel = nelem (); |
4634 | 582 |
5275 | 583 for (octave_idx_type i = 0; i < nel; i++) |
4634 | 584 { |
585 double val = elem (i); | |
586 if (xisinf (val) || xisnan (val)) | |
587 return true; | |
588 } | |
589 | |
590 return false; | |
591 } | |
592 | |
593 bool | |
5943 | 594 NDArray::any_element_not_one_or_zero (void) const |
595 { | |
596 octave_idx_type nel = nelem (); | |
597 | |
598 for (octave_idx_type i = 0; i < nel; i++) | |
599 { | |
600 double val = elem (i); | |
601 if (val != 0 && val != 1) | |
602 return true; | |
603 } | |
604 | |
605 return false; | |
606 } | |
607 | |
608 bool | |
6989 | 609 NDArray::all_elements_are_zero (void) const |
610 { | |
611 octave_idx_type nel = nelem (); | |
612 | |
613 for (octave_idx_type i = 0; i < nel; i++) | |
614 if (elem (i) != 0) | |
615 return false; | |
616 | |
617 return true; | |
618 } | |
619 | |
620 bool | |
4634 | 621 NDArray::all_elements_are_int_or_inf_or_nan (void) const |
622 { | |
5275 | 623 octave_idx_type nel = nelem (); |
4634 | 624 |
5275 | 625 for (octave_idx_type i = 0; i < nel; i++) |
4634 | 626 { |
627 double val = elem (i); | |
628 if (xisnan (val) || D_NINT (val) == val) | |
629 continue; | |
630 else | |
631 return false; | |
632 } | |
633 | |
634 return true; | |
635 } | |
636 | |
637 // Return nonzero if any element of M is not an integer. Also extract | |
638 // the largest and smallest values and return them in MAX_VAL and MIN_VAL. | |
639 | |
640 bool | |
641 NDArray::all_integers (double& max_val, double& min_val) const | |
642 { | |
5275 | 643 octave_idx_type nel = nelem (); |
4634 | 644 |
645 if (nel > 0) | |
646 { | |
647 max_val = elem (0); | |
648 min_val = elem (0); | |
649 } | |
650 else | |
651 return false; | |
652 | |
5275 | 653 for (octave_idx_type i = 0; i < nel; i++) |
4634 | 654 { |
655 double val = elem (i); | |
656 | |
657 if (val > max_val) | |
658 max_val = val; | |
659 | |
660 if (val < min_val) | |
661 min_val = val; | |
662 | |
663 if (D_NINT (val) != val) | |
664 return false; | |
665 } | |
666 | |
667 return true; | |
668 } | |
669 | |
670 bool | |
671 NDArray::too_large_for_float (void) const | |
672 { | |
5275 | 673 octave_idx_type nel = nelem (); |
4634 | 674 |
5275 | 675 for (octave_idx_type i = 0; i < nel; i++) |
4634 | 676 { |
677 double val = elem (i); | |
678 | |
5389 | 679 if (! (xisnan (val) || xisinf (val)) |
5387 | 680 && fabs (val) > FLT_MAX) |
4634 | 681 return true; |
682 } | |
683 | |
684 return false; | |
685 } | |
686 | |
5775 | 687 // FIXME -- this is not quite the right thing. |
4513 | 688 |
4556 | 689 boolNDArray |
4513 | 690 NDArray::all (int dim) const |
691 { | |
8743
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
692 return do_mx_red_op<boolNDArray> (*this, dim, mx_inline_all); |
4513 | 693 } |
694 | |
4556 | 695 boolNDArray |
4513 | 696 NDArray::any (int dim) const |
697 { | |
8743
1bd918cfb6e2
reimplement any & all using the new reduction code
Jaroslav Hajek <highegg@gmail.com>
parents:
8736
diff
changeset
|
698 return do_mx_red_op<boolNDArray> (*this, dim, mx_inline_any); |
4569 | 699 } |
700 | |
4584 | 701 NDArray |
4569 | 702 NDArray::cumprod (int dim) const |
703 { | |
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
704 return do_mx_cum_op<NDArray> (*this, dim, mx_inline_cumprod); |
4569 | 705 } |
706 | |
4584 | 707 NDArray |
4569 | 708 NDArray::cumsum (int dim) const |
709 { | |
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
710 return do_mx_cum_op<NDArray> (*this, dim, mx_inline_cumsum); |
4513 | 711 } |
712 | |
4569 | 713 NDArray |
714 NDArray::prod (int dim) const | |
715 { | |
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
716 return do_mx_red_op<NDArray> (*this, dim, mx_inline_prod); |
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
717 } |
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
718 |
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
719 NDArray |
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
720 NDArray::sum (int dim) const |
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
721 { |
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
722 return do_mx_red_op<NDArray> (*this, dim, mx_inline_sum); |
4569 | 723 } |
724 | |
725 NDArray | |
726 NDArray::sumsq (int dim) const | |
727 { | |
8736
53b4fdeacc2e
improve reduction functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8650
diff
changeset
|
728 return do_mx_red_op<NDArray> (*this, dim, mx_inline_sumsq); |
4569 | 729 } |
730 | |
4844 | 731 NDArray |
732 NDArray::max (int dim) const | |
733 { | |
8751
9f7ce4bf7650
optimize min/max functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8743
diff
changeset
|
734 return do_mx_minmax_op<NDArray> (*this, dim, mx_inline_max); |
4844 | 735 } |
736 | |
737 NDArray | |
5275 | 738 NDArray::max (ArrayN<octave_idx_type>& idx_arg, int dim) const |
4844 | 739 { |
8751
9f7ce4bf7650
optimize min/max functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8743
diff
changeset
|
740 return do_mx_minmax_op<NDArray> (*this, idx_arg, dim, mx_inline_max); |
4844 | 741 } |
742 | |
743 NDArray | |
744 NDArray::min (int dim) const | |
745 { | |
8751
9f7ce4bf7650
optimize min/max functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8743
diff
changeset
|
746 return do_mx_minmax_op<NDArray> (*this, dim, mx_inline_min); |
4844 | 747 } |
748 | |
749 NDArray | |
5275 | 750 NDArray::min (ArrayN<octave_idx_type>& idx_arg, int dim) const |
4844 | 751 { |
8751
9f7ce4bf7650
optimize min/max functions
Jaroslav Hajek <highegg@gmail.com>
parents:
8743
diff
changeset
|
752 return do_mx_minmax_op<NDArray> (*this, idx_arg, dim, mx_inline_min); |
4844 | 753 } |
754 | |
4915 | 755 NDArray |
5275 | 756 NDArray::concat (const NDArray& rb, const Array<octave_idx_type>& ra_idx) |
4758 | 757 { |
5073 | 758 if (rb.numel () > 0) |
759 insert (rb, ra_idx); | |
760 return *this; | |
761 } | |
762 | |
763 ComplexNDArray | |
5275 | 764 NDArray::concat (const ComplexNDArray& rb, const Array<octave_idx_type>& ra_idx) |
5073 | 765 { |
766 ComplexNDArray retval (*this); | |
4940 | 767 if (rb.numel () > 0) |
4915 | 768 retval.insert (rb, ra_idx); |
769 return retval; | |
4758 | 770 } |
771 | |
5073 | 772 charNDArray |
5275 | 773 NDArray::concat (const charNDArray& rb, const Array<octave_idx_type>& ra_idx) |
5073 | 774 { |
775 charNDArray retval (dims ()); | |
5275 | 776 octave_idx_type nel = numel (); |
5073 | 777 |
5275 | 778 for (octave_idx_type i = 0; i < nel; i++) |
5073 | 779 { |
780 double d = elem (i); | |
781 | |
782 if (xisnan (d)) | |
783 { | |
784 (*current_liboctave_error_handler) | |
785 ("invalid conversion from NaN to character"); | |
786 return retval; | |
787 } | |
788 else | |
789 { | |
5275 | 790 octave_idx_type ival = NINTbig (d); |
5073 | 791 |
792 if (ival < 0 || ival > UCHAR_MAX) | |
5775 | 793 // FIXME -- is there something |
5073 | 794 // better we could do? Should we warn the user? |
795 ival = 0; | |
796 | |
797 retval.elem (i) = static_cast<char>(ival); | |
798 } | |
799 } | |
800 | |
801 if (rb.numel () == 0) | |
802 return retval; | |
803 | |
804 retval.insert (rb, ra_idx); | |
805 return retval; | |
806 } | |
807 | |
4634 | 808 NDArray |
809 real (const ComplexNDArray& a) | |
810 { | |
8650
a1ae2aae903e
abs,real,imag,conj: use code from mx-inlines rather than the generic map
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
811 return NDArray (mx_inline_real_dup (a.data (), a.length ()), |
a1ae2aae903e
abs,real,imag,conj: use code from mx-inlines rather than the generic map
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
812 a.dims ()); |
4634 | 813 } |
814 | |
815 NDArray | |
816 imag (const ComplexNDArray& a) | |
817 { | |
8650
a1ae2aae903e
abs,real,imag,conj: use code from mx-inlines rather than the generic map
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
818 return NDArray (mx_inline_imag_dup (a.data (), a.length ()), |
a1ae2aae903e
abs,real,imag,conj: use code from mx-inlines rather than the generic map
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
819 a.dims ()); |
4634 | 820 } |
821 | |
4915 | 822 NDArray& |
5275 | 823 NDArray::insert (const NDArray& a, octave_idx_type r, octave_idx_type c) |
4915 | 824 { |
825 Array<double>::insert (a, r, c); | |
826 return *this; | |
827 } | |
828 | |
829 NDArray& | |
5275 | 830 NDArray::insert (const NDArray& a, const Array<octave_idx_type>& ra_idx) |
4915 | 831 { |
832 Array<double>::insert (a, ra_idx); | |
833 return *this; | |
834 } | |
835 | |
4634 | 836 NDArray |
4569 | 837 NDArray::abs (void) const |
838 { | |
8650
a1ae2aae903e
abs,real,imag,conj: use code from mx-inlines rather than the generic map
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
839 return NDArray (mx_inline_fabs_dup (data (), length ()), |
a1ae2aae903e
abs,real,imag,conj: use code from mx-inlines rather than the generic map
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
840 dims ()); |
4569 | 841 } |
842 | |
4532 | 843 Matrix |
844 NDArray::matrix_value (void) const | |
845 { | |
846 Matrix retval; | |
847 | |
848 int nd = ndims (); | |
849 | |
850 switch (nd) | |
851 { | |
852 case 1: | |
853 retval = Matrix (Array2<double> (*this, dimensions(0), 1)); | |
854 break; | |
855 | |
856 case 2: | |
857 retval = Matrix (Array2<double> (*this, dimensions(0), dimensions(1))); | |
858 break; | |
859 | |
860 default: | |
861 (*current_liboctave_error_handler) | |
4770 | 862 ("invalid conversion of NDArray to Matrix"); |
4532 | 863 break; |
864 } | |
865 | |
866 return retval; | |
867 } | |
868 | |
869 void | |
5275 | 870 NDArray::increment_index (Array<octave_idx_type>& ra_idx, |
4532 | 871 const dim_vector& dimensions, |
872 int start_dimension) | |
873 { | |
874 ::increment_index (ra_idx, dimensions, start_dimension); | |
875 } | |
876 | |
5275 | 877 octave_idx_type |
878 NDArray::compute_index (Array<octave_idx_type>& ra_idx, | |
4556 | 879 const dim_vector& dimensions) |
880 { | |
881 return ::compute_index (ra_idx, dimensions); | |
882 } | |
883 | |
7503
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
884 NDArray |
7620
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7600
diff
changeset
|
885 NDArray::diag (octave_idx_type k) const |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7600
diff
changeset
|
886 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7600
diff
changeset
|
887 return MArrayN<double>::diag (k); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7600
diff
changeset
|
888 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7600
diff
changeset
|
889 |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7600
diff
changeset
|
890 NDArray |
7503
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
891 NDArray::map (dmapper fcn) const |
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
892 { |
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
893 return MArrayN<double>::map<double> (func_ptr (fcn)); |
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
894 } |
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
895 |
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
896 ComplexNDArray |
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
897 NDArray::map (cmapper fcn) const |
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
898 { |
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
899 return MArrayN<double>::map<Complex> (func_ptr (fcn)); |
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
900 } |
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
901 |
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
902 boolNDArray |
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
903 NDArray::map (bmapper fcn) const |
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
904 { |
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
905 return MArrayN<double>::map<bool> (func_ptr (fcn)); |
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
906 } |
8c32f95c2639
convert mapper functions to new format
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
907 |
4687 | 908 // This contains no information on the array structure !!! |
909 std::ostream& | |
910 operator << (std::ostream& os, const NDArray& a) | |
911 { | |
5275 | 912 octave_idx_type nel = a.nelem (); |
4687 | 913 |
5275 | 914 for (octave_idx_type i = 0; i < nel; i++) |
4687 | 915 { |
916 os << " "; | |
917 octave_write_double (os, a.elem (i)); | |
918 os << "\n"; | |
919 } | |
920 return os; | |
921 } | |
922 | |
923 std::istream& | |
924 operator >> (std::istream& is, NDArray& a) | |
925 { | |
5275 | 926 octave_idx_type nel = a.nelem (); |
4687 | 927 |
928 if (nel < 1 ) | |
929 is.clear (std::ios::badbit); | |
930 else | |
931 { | |
932 double tmp; | |
5275 | 933 for (octave_idx_type i = 0; i < nel; i++) |
4687 | 934 { |
935 tmp = octave_read_double (is); | |
936 if (is) | |
937 a.elem (i) = tmp; | |
938 else | |
939 goto done; | |
940 } | |
941 } | |
942 | |
943 done: | |
944 | |
945 return is; | |
946 } | |
947 | |
5775 | 948 // FIXME -- it would be nice to share code among the min/max |
4844 | 949 // functions below. |
950 | |
951 #define EMPTY_RETURN_CHECK(T) \ | |
952 if (nel == 0) \ | |
953 return T (dv); | |
954 | |
955 NDArray | |
956 min (double d, const NDArray& m) | |
957 { | |
958 dim_vector dv = m.dims (); | |
5275 | 959 octave_idx_type nel = dv.numel (); |
4844 | 960 |
961 EMPTY_RETURN_CHECK (NDArray); | |
962 | |
963 NDArray result (dv); | |
964 | |
5275 | 965 for (octave_idx_type i = 0; i < nel; i++) |
4844 | 966 { |
967 OCTAVE_QUIT; | |
968 result (i) = xmin (d, m (i)); | |
969 } | |
970 | |
971 return result; | |
972 } | |
973 | |
974 NDArray | |
975 min (const NDArray& m, double d) | |
976 { | |
977 dim_vector dv = m.dims (); | |
5275 | 978 octave_idx_type nel = dv.numel (); |
4844 | 979 |
980 EMPTY_RETURN_CHECK (NDArray); | |
981 | |
982 NDArray result (dv); | |
983 | |
5275 | 984 for (octave_idx_type i = 0; i < nel; i++) |
4844 | 985 { |
986 OCTAVE_QUIT; | |
987 result (i) = xmin (d, m (i)); | |
988 } | |
989 | |
990 return result; | |
991 } | |
992 | |
993 NDArray | |
994 min (const NDArray& a, const NDArray& b) | |
995 { | |
996 dim_vector dv = a.dims (); | |
5275 | 997 octave_idx_type nel = dv.numel (); |
4844 | 998 |
999 if (dv != b.dims ()) | |
1000 { | |
1001 (*current_liboctave_error_handler) | |
1002 ("two-arg min expecting args of same size"); | |
1003 return NDArray (); | |
1004 } | |
1005 | |
1006 EMPTY_RETURN_CHECK (NDArray); | |
1007 | |
1008 NDArray result (dv); | |
1009 | |
5275 | 1010 for (octave_idx_type i = 0; i < nel; i++) |
4844 | 1011 { |
1012 OCTAVE_QUIT; | |
1013 result (i) = xmin (a (i), b (i)); | |
1014 } | |
1015 | |
1016 return result; | |
1017 } | |
1018 | |
1019 NDArray | |
1020 max (double d, const NDArray& m) | |
1021 { | |
1022 dim_vector dv = m.dims (); | |
5275 | 1023 octave_idx_type nel = dv.numel (); |
4844 | 1024 |
1025 EMPTY_RETURN_CHECK (NDArray); | |
1026 | |
1027 NDArray result (dv); | |
1028 | |
5275 | 1029 for (octave_idx_type i = 0; i < nel; i++) |
4844 | 1030 { |
1031 OCTAVE_QUIT; | |
1032 result (i) = xmax (d, m (i)); | |
1033 } | |
1034 | |
1035 return result; | |
1036 } | |
1037 | |
1038 NDArray | |
1039 max (const NDArray& m, double d) | |
1040 { | |
1041 dim_vector dv = m.dims (); | |
5275 | 1042 octave_idx_type nel = dv.numel (); |
4844 | 1043 |
1044 EMPTY_RETURN_CHECK (NDArray); | |
1045 | |
1046 NDArray result (dv); | |
1047 | |
5275 | 1048 for (octave_idx_type i = 0; i < nel; i++) |
4844 | 1049 { |
1050 OCTAVE_QUIT; | |
1051 result (i) = xmax (d, m (i)); | |
1052 } | |
1053 | |
1054 return result; | |
1055 } | |
1056 | |
1057 NDArray | |
1058 max (const NDArray& a, const NDArray& b) | |
1059 { | |
1060 dim_vector dv = a.dims (); | |
5275 | 1061 octave_idx_type nel = dv.numel (); |
4844 | 1062 |
1063 if (dv != b.dims ()) | |
1064 { | |
1065 (*current_liboctave_error_handler) | |
1066 ("two-arg max expecting args of same size"); | |
1067 return NDArray (); | |
1068 } | |
1069 | |
1070 EMPTY_RETURN_CHECK (NDArray); | |
1071 | |
1072 NDArray result (dv); | |
1073 | |
5275 | 1074 for (octave_idx_type i = 0; i < nel; i++) |
4844 | 1075 { |
1076 OCTAVE_QUIT; | |
1077 result (i) = xmax (a (i), b (i)); | |
1078 } | |
1079 | |
1080 return result; | |
1081 } | |
1082 | |
4543 | 1083 NDS_CMP_OPS(NDArray, , double, ) |
1084 NDS_BOOL_OPS(NDArray, double, 0.0) | |
1085 | |
1086 SND_CMP_OPS(double, , NDArray, ) | |
1087 SND_BOOL_OPS(double, NDArray, 0.0) | |
1088 | |
1089 NDND_CMP_OPS(NDArray, , NDArray, ) | |
1090 NDND_BOOL_OPS(NDArray, NDArray, 0.0) | |
1091 | |
4513 | 1092 /* |
1093 ;;; Local Variables: *** | |
1094 ;;; mode: C++ *** | |
1095 ;;; End: *** | |
1096 */ |