Mercurial > hg > octave-lyh

comparison liboctave/dNDArray.cc @ 5275:23b37da9fd5b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-04-08 16:07:35 by jwe]
author jwe
date Fri, 08 Apr 2005 16:07:37 +0000
parents 57077d0ddc8e
children 4c8a2e4e0717
comparison
equal deleted inserted replaced
5274:eae7b40388e9 5275:23b37da9fd5b
46 dim_vector dv = dims (); 46 dim_vector dv = dims ();
47 47
48 if (dim > dv.length () || dim < 0) 48 if (dim > dv.length () || dim < 0)
49 return ComplexNDArray (); 49 return ComplexNDArray ();
50 50
51 int stride = 1; 51 octave_idx_type stride = 1;
52 int n = dv(dim); 52 octave_idx_type n = dv(dim);
53 53
54 for (int i = 0; i < dim; i++) 54 for (int i = 0; i < dim; i++)
55 stride *= dv(i); 55 stride *= dv(i);
56 56
57 int howmany = numel () / dv (dim); 57 octave_idx_type howmany = numel () / dv (dim);
58 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany)); 58 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
59 int nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride); 59 octave_idx_type nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride);
60 int dist = (stride == 1 ? n : 1); 60 octave_idx_type dist = (stride == 1 ? n : 1);
61 61
62 const double *in (fortran_vec ()); 62 const double *in (fortran_vec ());
63 ComplexNDArray retval (dv); 63 ComplexNDArray retval (dv);
64 Complex *out (retval.fortran_vec ()); 64 Complex *out (retval.fortran_vec ());
65 65
66 // Need to be careful here about the distance between fft's 66 // Need to be careful here about the distance between fft's
67 for (int k = 0; k < nloop; k++) 67 for (octave_idx_type k = 0; k < nloop; k++)
68 octave_fftw::fft (in + k * stride * n, out + k * stride * n, 68 octave_fftw::fft (in + k * stride * n, out + k * stride * n,
69 n, howmany, stride, dist); 69 n, howmany, stride, dist);
70 70
71 return retval; 71 return retval;
72 } 72 }
77 dim_vector dv = dims (); 77 dim_vector dv = dims ();
78 78
79 if (dim > dv.length () || dim < 0) 79 if (dim > dv.length () || dim < 0)
80 return ComplexNDArray (); 80 return ComplexNDArray ();
81 81
82 int stride = 1; 82 octave_idx_type stride = 1;
83 int n = dv(dim); 83 octave_idx_type n = dv(dim);
84 84
85 for (int i = 0; i < dim; i++) 85 for (int i = 0; i < dim; i++)
86 stride *= dv(i); 86 stride *= dv(i);
87 87
88 int howmany = numel () / dv (dim); 88 octave_idx_type howmany = numel () / dv (dim);
89 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany)); 89 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
90 int nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride); 90 octave_idx_type nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride);
91 int dist = (stride == 1 ? n : 1); 91 octave_idx_type dist = (stride == 1 ? n : 1);
92 92
93 ComplexNDArray retval (*this); 93 ComplexNDArray retval (*this);
94 Complex *out (retval.fortran_vec ()); 94 Complex *out (retval.fortran_vec ());
95 95
96 // Need to be careful here about the distance between fft's 96 // Need to be careful here about the distance between fft's
97 for (int k = 0; k < nloop; k++) 97 for (octave_idx_type k = 0; k < nloop; k++)
98 octave_fftw::ifft (out + k * stride * n, out + k * stride * n, 98 octave_fftw::ifft (out + k * stride * n, out + k * stride * n,
99 n, howmany, stride, dist); 99 n, howmany, stride, dist);
100 100
101 return retval; 101 return retval;
102 } 102 }
110 110
111 dim_vector dv2(dv(0), dv(1)); 111 dim_vector dv2(dv(0), dv(1));
112 const double *in = fortran_vec (); 112 const double *in = fortran_vec ();
113 ComplexNDArray retval (dv); 113 ComplexNDArray retval (dv);
114 Complex *out = retval.fortran_vec (); 114 Complex *out = retval.fortran_vec ();
115 int howmany = numel() / dv(0) / dv(1); 115 octave_idx_type howmany = numel() / dv(0) / dv(1);
116 int dist = dv(0) * dv(1); 116 octave_idx_type dist = dv(0) * dv(1);
117 117
118 for (int i=0; i < howmany; i++) 118 for (octave_idx_type i=0; i < howmany; i++)
119 octave_fftw::fftNd (in + i*dist, out + i*dist, 2, dv2); 119 octave_fftw::fftNd (in + i*dist, out + i*dist, 2, dv2);
120 120
121 return retval; 121 return retval;
122 } 122 }
123 123
129 return ComplexNDArray (); 129 return ComplexNDArray ();
130 130
131 dim_vector dv2(dv(0), dv(1)); 131 dim_vector dv2(dv(0), dv(1));
132 ComplexNDArray retval (*this); 132 ComplexNDArray retval (*this);
133 Complex *out = retval.fortran_vec (); 133 Complex *out = retval.fortran_vec ();
134 int howmany = numel() / dv(0) / dv(1); 134 octave_idx_type howmany = numel() / dv(0) / dv(1);
135 int dist = dv(0) * dv(1); 135 octave_idx_type dist = dv(0) * dv(1);
136 136
137 for (int i=0; i < howmany; i++) 137 for (octave_idx_type i=0; i < howmany; i++)
138 octave_fftw::ifftNd (out + i*dist, out + i*dist, 2, dv2); 138 octave_fftw::ifftNd (out + i*dist, out + i*dist, 2, dv2);
139 139
140 return retval; 140 return retval;
141 } 141 }
142 142
179 // double complex arguments. They have been modified by adding an 179 // double complex arguments. They have been modified by adding an
180 // implicit double precision (a-h,o-z) statement at the beginning of 180 // implicit double precision (a-h,o-z) statement at the beginning of
181 // each subroutine. 181 // each subroutine.
182 182
183 F77_RET_T 183 F77_RET_T
184 F77_FUNC (cffti, CFFTI) (const int&, Complex*); 184 F77_FUNC (cffti, CFFTI) (const octave_idx_type&, Complex*);
185 185
186 F77_RET_T 186 F77_RET_T
187 F77_FUNC (cfftf, CFFTF) (const int&, Complex*, Complex*); 187 F77_FUNC (cfftf, CFFTF) (const octave_idx_type&, Complex*, Complex*);
188 188
189 F77_RET_T 189 F77_RET_T
190 F77_FUNC (cfftb, CFFTB) (const int&, Complex*, Complex*); 190 F77_FUNC (cfftb, CFFTB) (const octave_idx_type&, Complex*, Complex*);
191 } 191 }
192 192
193 ComplexNDArray 193 ComplexNDArray
194 NDArray::fourier (int dim) const 194 NDArray::fourier (int dim) const
195 { 195 {
197 197
198 if (dim > dv.length () || dim < 0) 198 if (dim > dv.length () || dim < 0)
199 return ComplexNDArray (); 199 return ComplexNDArray ();
200 200
201 ComplexNDArray retval (dv); 201 ComplexNDArray retval (dv);
202 int npts = dv(dim); 202 octave_idx_type npts = dv(dim);
203 int nn = 4*npts+15; 203 octave_idx_type nn = 4*npts+15;
204 Array<Complex> wsave (nn); 204 Array<Complex> wsave (nn);
205 Complex *pwsave = wsave.fortran_vec (); 205 Complex *pwsave = wsave.fortran_vec ();
206 206
207 OCTAVE_LOCAL_BUFFER (Complex, tmp, npts); 207 OCTAVE_LOCAL_BUFFER (Complex, tmp, npts);
208 208
209 int stride = 1; 209 octave_idx_type stride = 1;
210 210
211 for (int i = 0; i < dim; i++) 211 for (int i = 0; i < dim; i++)
212 stride *= dv(i); 212 stride *= dv(i);
213 213
214 int howmany = numel () / npts; 214 octave_idx_type howmany = numel () / npts;
215 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany)); 215 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
216 int nloop = (stride == 1 ? 1 : numel () / npts / stride); 216 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride);
217 int dist = (stride == 1 ? npts : 1); 217 octave_idx_type dist = (stride == 1 ? npts : 1);
218 218
219 F77_FUNC (cffti, CFFTI) (npts, pwsave); 219 F77_FUNC (cffti, CFFTI) (npts, pwsave);
220 220
221 for (int k = 0; k < nloop; k++) 221 for (octave_idx_type k = 0; k < nloop; k++)
222 { 222 {
223 for (int j = 0; j < howmany; j++) 223 for (octave_idx_type j = 0; j < howmany; j++)
224 { 224 {
225 OCTAVE_QUIT; 225 OCTAVE_QUIT;
226 226
227 for (int i = 0; i < npts; i++) 227 for (octave_idx_type i = 0; i < npts; i++)
228 tmp[i] = elem((i + k*npts)*stride + j*dist); 228 tmp[i] = elem((i + k*npts)*stride + j*dist);
229 229
230 F77_FUNC (cfftf, CFFTF) (npts, tmp, pwsave); 230 F77_FUNC (cfftf, CFFTF) (npts, tmp, pwsave);
231 231
232 for (int i = 0; i < npts; i++) 232 for (octave_idx_type i = 0; i < npts; i++)
233 retval ((i + k*npts)*stride + j*dist) = tmp[i]; 233 retval ((i + k*npts)*stride + j*dist) = tmp[i];
234 } 234 }
235 } 235 }
236 236
237 return retval; 237 return retval;
244 244
245 if (dim > dv.length () || dim < 0) 245 if (dim > dv.length () || dim < 0)
246 return ComplexNDArray (); 246 return ComplexNDArray ();
247 247
248 ComplexNDArray retval (dv); 248 ComplexNDArray retval (dv);
249 int npts = dv(dim); 249 octave_idx_type npts = dv(dim);
250 int nn = 4*npts+15; 250 octave_idx_type nn = 4*npts+15;
251 Array<Complex> wsave (nn); 251 Array<Complex> wsave (nn);
252 Complex *pwsave = wsave.fortran_vec (); 252 Complex *pwsave = wsave.fortran_vec ();
253 253
254 OCTAVE_LOCAL_BUFFER (Complex, tmp, npts); 254 OCTAVE_LOCAL_BUFFER (Complex, tmp, npts);
255 255
256 int stride = 1; 256 octave_idx_type stride = 1;
257 257
258 for (int i = 0; i < dim; i++) 258 for (int i = 0; i < dim; i++)
259 stride *= dv(i); 259 stride *= dv(i);
260 260
261 int howmany = numel () / npts; 261 octave_idx_type howmany = numel () / npts;
262 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany)); 262 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
263 int nloop = (stride == 1 ? 1 : numel () / npts / stride); 263 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride);
264 int dist = (stride == 1 ? npts : 1); 264 octave_idx_type dist = (stride == 1 ? npts : 1);
265 265
266 F77_FUNC (cffti, CFFTI) (npts, pwsave); 266 F77_FUNC (cffti, CFFTI) (npts, pwsave);
267 267
268 for (int k = 0; k < nloop; k++) 268 for (octave_idx_type k = 0; k < nloop; k++)
269 { 269 {
270 for (int j = 0; j < howmany; j++) 270 for (octave_idx_type j = 0; j < howmany; j++)
271 { 271 {
272 OCTAVE_QUIT; 272 OCTAVE_QUIT;
273 273
274 for (int i = 0; i < npts; i++) 274 for (octave_idx_type i = 0; i < npts; i++)
275 tmp[i] = elem((i + k*npts)*stride + j*dist); 275 tmp[i] = elem((i + k*npts)*stride + j*dist);
276 276
277 F77_FUNC (cfftb, CFFTB) (npts, tmp, pwsave); 277 F77_FUNC (cfftb, CFFTB) (npts, tmp, pwsave);
278 278
279 for (int i = 0; i < npts; i++) 279 for (octave_idx_type i = 0; i < npts; i++)
280 retval ((i + k*npts)*stride + j*dist) = tmp[i] / 280 retval ((i + k*npts)*stride + j*dist) = tmp[i] /
281 static_cast<double> (npts); 281 static_cast<double> (npts);
282 } 282 }
283 } 283 }
284 284
290 { 290 {
291 dim_vector dv = dims(); 291 dim_vector dv = dims();
292 dim_vector dv2 (dv(0), dv(1)); 292 dim_vector dv2 (dv(0), dv(1));
293 int rank = 2; 293 int rank = 2;
294 ComplexNDArray retval (*this); 294 ComplexNDArray retval (*this);
295 int stride = 1; 295 octave_idx_type stride = 1;
296 296
297 for (int i = 0; i < rank; i++) 297 for (int i = 0; i < rank; i++)
298 { 298 {
299 int npts = dv2(i); 299 octave_idx_type npts = dv2(i);
300 int nn = 4*npts+15; 300 octave_idx_type nn = 4*npts+15;
301 Array<Complex> wsave (nn); 301 Array<Complex> wsave (nn);
302 Complex *pwsave = wsave.fortran_vec (); 302 Complex *pwsave = wsave.fortran_vec ();
303 Array<Complex> row (npts); 303 Array<Complex> row (npts);
304 Complex *prow = row.fortran_vec (); 304 Complex *prow = row.fortran_vec ();
305 305
306 int howmany = numel () / npts; 306 octave_idx_type howmany = numel () / npts;
307 howmany = (stride == 1 ? howmany : 307 howmany = (stride == 1 ? howmany :
308 (howmany > stride ? stride : howmany)); 308 (howmany > stride ? stride : howmany));
309 int nloop = (stride == 1 ? 1 : numel () / npts / stride); 309 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride);
310 int dist = (stride == 1 ? npts : 1); 310 octave_idx_type dist = (stride == 1 ? npts : 1);
311 311
312 F77_FUNC (cffti, CFFTI) (npts, pwsave); 312 F77_FUNC (cffti, CFFTI) (npts, pwsave);
313 313
314 for (int k = 0; k < nloop; k++) 314 for (octave_idx_type k = 0; k < nloop; k++)
315 { 315 {
316 for (int j = 0; j < howmany; j++) 316 for (octave_idx_type j = 0; j < howmany; j++)
317 { 317 {
318 OCTAVE_QUIT; 318 OCTAVE_QUIT;
319 319
320 for (int l = 0; l < npts; l++) 320 for (octave_idx_type l = 0; l < npts; l++)
321 prow[l] = retval ((l + k*npts)*stride + j*dist); 321 prow[l] = retval ((l + k*npts)*stride + j*dist);
322 322
323 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave); 323 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
324 324
325 for (int l = 0; l < npts; l++) 325 for (octave_idx_type l = 0; l < npts; l++)
326 retval ((l + k*npts)*stride + j*dist) = prow[l]; 326 retval ((l + k*npts)*stride + j*dist) = prow[l];
327 } 327 }
328 } 328 }
329 329
330 stride *= dv2(i); 330 stride *= dv2(i);
338 { 338 {
339 dim_vector dv = dims(); 339 dim_vector dv = dims();
340 dim_vector dv2 (dv(0), dv(1)); 340 dim_vector dv2 (dv(0), dv(1));
341 int rank = 2; 341 int rank = 2;
342 ComplexNDArray retval (*this); 342 ComplexNDArray retval (*this);
343 int stride = 1; 343 octave_idx_type stride = 1;
344 344
345 for (int i = 0; i < rank; i++) 345 for (int i = 0; i < rank; i++)
346 { 346 {
347 int npts = dv2(i); 347 octave_idx_type npts = dv2(i);
348 int nn = 4*npts+15; 348 octave_idx_type nn = 4*npts+15;
349 Array<Complex> wsave (nn); 349 Array<Complex> wsave (nn);
350 Complex *pwsave = wsave.fortran_vec (); 350 Complex *pwsave = wsave.fortran_vec ();
351 Array<Complex> row (npts); 351 Array<Complex> row (npts);
352 Complex *prow = row.fortran_vec (); 352 Complex *prow = row.fortran_vec ();
353 353
354 int howmany = numel () / npts; 354 octave_idx_type howmany = numel () / npts;
355 howmany = (stride == 1 ? howmany : 355 howmany = (stride == 1 ? howmany :
356 (howmany > stride ? stride : howmany)); 356 (howmany > stride ? stride : howmany));
357 int nloop = (stride == 1 ? 1 : numel () / npts / stride); 357 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride);
358 int dist = (stride == 1 ? npts : 1); 358 octave_idx_type dist = (stride == 1 ? npts : 1);
359 359
360 F77_FUNC (cffti, CFFTI) (npts, pwsave); 360 F77_FUNC (cffti, CFFTI) (npts, pwsave);
361 361
362 for (int k = 0; k < nloop; k++) 362 for (octave_idx_type k = 0; k < nloop; k++)
363 { 363 {
364 for (int j = 0; j < howmany; j++) 364 for (octave_idx_type j = 0; j < howmany; j++)
365 { 365 {
366 OCTAVE_QUIT; 366 OCTAVE_QUIT;
367 367
368 for (int l = 0; l < npts; l++) 368 for (octave_idx_type l = 0; l < npts; l++)
369 prow[l] = retval ((l + k*npts)*stride + j*dist); 369 prow[l] = retval ((l + k*npts)*stride + j*dist);
370 370
371 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave); 371 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
372 372
373 for (int l = 0; l < npts; l++) 373 for (octave_idx_type l = 0; l < npts; l++)
374 retval ((l + k*npts)*stride + j*dist) = prow[l] / 374 retval ((l + k*npts)*stride + j*dist) = prow[l] /
375 static_cast<double> (npts); 375 static_cast<double> (npts);
376 } 376 }
377 } 377 }
378 378
386 NDArray::fourierNd (void) const 386 NDArray::fourierNd (void) const
387 { 387 {
388 dim_vector dv = dims (); 388 dim_vector dv = dims ();
389 int rank = dv.length (); 389 int rank = dv.length ();
390 ComplexNDArray retval (*this); 390 ComplexNDArray retval (*this);
391 int stride = 1; 391 octave_idx_type stride = 1;
392 392
393 for (int i = 0; i < rank; i++) 393 for (int i = 0; i < rank; i++)
394 { 394 {
395 int npts = dv(i); 395 octave_idx_type npts = dv(i);
396 int nn = 4*npts+15; 396 octave_idx_type nn = 4*npts+15;
397 Array<Complex> wsave (nn); 397 Array<Complex> wsave (nn);
398 Complex *pwsave = wsave.fortran_vec (); 398 Complex *pwsave = wsave.fortran_vec ();
399 Array<Complex> row (npts); 399 Array<Complex> row (npts);
400 Complex *prow = row.fortran_vec (); 400 Complex *prow = row.fortran_vec ();
401 401
402 int howmany = numel () / npts; 402 octave_idx_type howmany = numel () / npts;
403 howmany = (stride == 1 ? howmany : 403 howmany = (stride == 1 ? howmany :
404 (howmany > stride ? stride : howmany)); 404 (howmany > stride ? stride : howmany));
405 int nloop = (stride == 1 ? 1 : numel () / npts / stride); 405 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride);
406 int dist = (stride == 1 ? npts : 1); 406 octave_idx_type dist = (stride == 1 ? npts : 1);
407 407
408 F77_FUNC (cffti, CFFTI) (npts, pwsave); 408 F77_FUNC (cffti, CFFTI) (npts, pwsave);
409 409
410 for (int k = 0; k < nloop; k++) 410 for (octave_idx_type k = 0; k < nloop; k++)
411 { 411 {
412 for (int j = 0; j < howmany; j++) 412 for (octave_idx_type j = 0; j < howmany; j++)
413 { 413 {
414 OCTAVE_QUIT; 414 OCTAVE_QUIT;
415 415
416 for (int l = 0; l < npts; l++) 416 for (octave_idx_type l = 0; l < npts; l++)
417 prow[l] = retval ((l + k*npts)*stride + j*dist); 417 prow[l] = retval ((l + k*npts)*stride + j*dist);
418 418
419 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave); 419 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
420 420
421 for (int l = 0; l < npts; l++) 421 for (octave_idx_type l = 0; l < npts; l++)
422 retval ((l + k*npts)*stride + j*dist) = prow[l]; 422 retval ((l + k*npts)*stride + j*dist) = prow[l];
423 } 423 }
424 } 424 }
425 425
426 stride *= dv(i); 426 stride *= dv(i);
433 NDArray::ifourierNd (void) const 433 NDArray::ifourierNd (void) const
434 { 434 {
435 dim_vector dv = dims (); 435 dim_vector dv = dims ();
436 int rank = dv.length (); 436 int rank = dv.length ();
437 ComplexNDArray retval (*this); 437 ComplexNDArray retval (*this);
438 int stride = 1; 438 octave_idx_type stride = 1;
439 439
440 for (int i = 0; i < rank; i++) 440 for (int i = 0; i < rank; i++)
441 { 441 {
442 int npts = dv(i); 442 octave_idx_type npts = dv(i);
443 int nn = 4*npts+15; 443 octave_idx_type nn = 4*npts+15;
444 Array<Complex> wsave (nn); 444 Array<Complex> wsave (nn);
445 Complex *pwsave = wsave.fortran_vec (); 445 Complex *pwsave = wsave.fortran_vec ();
446 Array<Complex> row (npts); 446 Array<Complex> row (npts);
447 Complex *prow = row.fortran_vec (); 447 Complex *prow = row.fortran_vec ();
448 448
449 int howmany = numel () / npts; 449 octave_idx_type howmany = numel () / npts;
450 howmany = (stride == 1 ? howmany : 450 howmany = (stride == 1 ? howmany :
451 (howmany > stride ? stride : howmany)); 451 (howmany > stride ? stride : howmany));
452 int nloop = (stride == 1 ? 1 : numel () / npts / stride); 452 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride);
453 int dist = (stride == 1 ? npts : 1); 453 octave_idx_type dist = (stride == 1 ? npts : 1);
454 454
455 F77_FUNC (cffti, CFFTI) (npts, pwsave); 455 F77_FUNC (cffti, CFFTI) (npts, pwsave);
456 456
457 for (int k = 0; k < nloop; k++) 457 for (octave_idx_type k = 0; k < nloop; k++)
458 { 458 {
459 for (int j = 0; j < howmany; j++) 459 for (octave_idx_type j = 0; j < howmany; j++)
460 { 460 {
461 OCTAVE_QUIT; 461 OCTAVE_QUIT;
462 462
463 for (int l = 0; l < npts; l++) 463 for (octave_idx_type l = 0; l < npts; l++)
464 prow[l] = retval ((l + k*npts)*stride + j*dist); 464 prow[l] = retval ((l + k*npts)*stride + j*dist);
465 465
466 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave); 466 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
467 467
468 for (int l = 0; l < npts; l++) 468 for (octave_idx_type l = 0; l < npts; l++)
469 retval ((l + k*npts)*stride + j*dist) = prow[l] / 469 retval ((l + k*npts)*stride + j*dist) = prow[l] /
470 static_cast<double> (npts); 470 static_cast<double> (npts);
471 } 471 }
472 } 472 }
473 473
484 boolNDArray 484 boolNDArray
485 NDArray::operator ! (void) const 485 NDArray::operator ! (void) const
486 { 486 {
487 boolNDArray b (dims ()); 487 boolNDArray b (dims ());
488 488
489 for (int i = 0; i < length (); i++) 489 for (octave_idx_type i = 0; i < length (); i++)
490 b.elem (i) = ! elem (i); 490 b.elem (i) = ! elem (i);
491 491
492 return b; 492 return b;
493 } 493 }
494 494
495 bool 495 bool
496 NDArray::any_element_is_negative (bool neg_zero) const 496 NDArray::any_element_is_negative (bool neg_zero) const
497 { 497 {
498 int nel = nelem (); 498 octave_idx_type nel = nelem ();
499 499
500 if (neg_zero) 500 if (neg_zero)
501 { 501 {
502 for (int i = 0; i < nel; i++) 502 for (octave_idx_type i = 0; i < nel; i++)
503 if (lo_ieee_signbit (elem (i))) 503 if (lo_ieee_signbit (elem (i)))
504 return true; 504 return true;
505 } 505 }
506 else 506 else
507 { 507 {
508 for (int i = 0; i < nel; i++) 508 for (octave_idx_type i = 0; i < nel; i++)
509 if (elem (i) < 0) 509 if (elem (i) < 0)
510 return true; 510 return true;
511 } 511 }
512 512
513 return false; 513 return false;
515 515
516 516
517 bool 517 bool
518 NDArray::any_element_is_inf_or_nan (void) const 518 NDArray::any_element_is_inf_or_nan (void) const
519 { 519 {
520 int nel = nelem (); 520 octave_idx_type nel = nelem ();
521 521
522 for (int i = 0; i < nel; i++) 522 for (octave_idx_type i = 0; i < nel; i++)
523 { 523 {
524 double val = elem (i); 524 double val = elem (i);
525 if (xisinf (val) || xisnan (val)) 525 if (xisinf (val) || xisnan (val))
526 return true; 526 return true;
527 } 527 }
530 } 530 }
531 531
532 bool 532 bool
533 NDArray::all_elements_are_int_or_inf_or_nan (void) const 533 NDArray::all_elements_are_int_or_inf_or_nan (void) const
534 { 534 {
535 int nel = nelem (); 535 octave_idx_type nel = nelem ();
536 536
537 for (int i = 0; i < nel; i++) 537 for (octave_idx_type i = 0; i < nel; i++)
538 { 538 {
539 double val = elem (i); 539 double val = elem (i);
540 if (xisnan (val) || D_NINT (val) == val) 540 if (xisnan (val) || D_NINT (val) == val)
541 continue; 541 continue;
542 else 542 else
550 // the largest and smallest values and return them in MAX_VAL and MIN_VAL. 550 // the largest and smallest values and return them in MAX_VAL and MIN_VAL.
551 551
552 bool 552 bool
553 NDArray::all_integers (double& max_val, double& min_val) const 553 NDArray::all_integers (double& max_val, double& min_val) const
554 { 554 {
555 int nel = nelem (); 555 octave_idx_type nel = nelem ();
556 556
557 if (nel > 0) 557 if (nel > 0)
558 { 558 {
559 max_val = elem (0); 559 max_val = elem (0);
560 min_val = elem (0); 560 min_val = elem (0);
561 } 561 }
562 else 562 else
563 return false; 563 return false;
564 564
565 for (int i = 0; i < nel; i++) 565 for (octave_idx_type i = 0; i < nel; i++)
566 { 566 {
567 double val = elem (i); 567 double val = elem (i);
568 568
569 if (val > max_val) 569 if (val > max_val)
570 max_val = val; 570 max_val = val;
580 } 580 }
581 581
582 bool 582 bool
583 NDArray::too_large_for_float (void) const 583 NDArray::too_large_for_float (void) const
584 { 584 {
585 int nel = nelem (); 585 octave_idx_type nel = nelem ();
586 586
587 for (int i = 0; i < nel; i++) 587 for (octave_idx_type i = 0; i < nel; i++)
588 { 588 {
589 double val = elem (i); 589 double val = elem (i);
590 590
591 if (val > FLT_MAX || val < FLT_MIN) 591 if (val > FLT_MAX || val < FLT_MIN)
592 return true; 592 return true;
642 } 642 }
643 643
644 NDArray 644 NDArray
645 NDArray::max (int dim) const 645 NDArray::max (int dim) const
646 { 646 {
647 ArrayN<int> dummy_idx; 647 ArrayN<octave_idx_type> dummy_idx;
648 return max (dummy_idx, dim); 648 return max (dummy_idx, dim);
649 } 649 }
650 650
651 NDArray 651 NDArray
652 NDArray::max (ArrayN<int>& idx_arg, int dim) const 652 NDArray::max (ArrayN<octave_idx_type>& idx_arg, int dim) const
653 { 653 {
654 dim_vector dv = dims (); 654 dim_vector dv = dims ();
655 dim_vector dr = dims (); 655 dim_vector dr = dims ();
656 656
657 if (dv.numel () == 0 || dim > dv.length () || dim < 0) 657 if (dv.numel () == 0 || dim > dv.length () || dim < 0)
660 dr(dim) = 1; 660 dr(dim) = 1;
661 661
662 NDArray result (dr); 662 NDArray result (dr);
663 idx_arg.resize (dr); 663 idx_arg.resize (dr);
664 664
665 int x_stride = 1; 665 octave_idx_type x_stride = 1;
666 int x_len = dv(dim); 666 octave_idx_type x_len = dv(dim);
667 for (int i = 0; i < dim; i++) 667 for (int i = 0; i < dim; i++)
668 x_stride *= dv(i); 668 x_stride *= dv(i);
669 669
670 for (int i = 0; i < dr.numel (); i++) 670 for (octave_idx_type i = 0; i < dr.numel (); i++)
671 { 671 {
672 int x_offset; 672 octave_idx_type x_offset;
673 if (x_stride == 1) 673 if (x_stride == 1)
674 x_offset = i * x_len; 674 x_offset = i * x_len;
675 else 675 else
676 { 676 {
677 int x_offset2 = 0; 677 octave_idx_type x_offset2 = 0;
678 x_offset = i; 678 x_offset = i;
679 while (x_offset >= x_stride) 679 while (x_offset >= x_stride)
680 { 680 {
681 x_offset -= x_stride; 681 x_offset -= x_stride;
682 x_offset2++; 682 x_offset2++;
683 } 683 }
684 x_offset += x_offset2 * x_stride * x_len; 684 x_offset += x_offset2 * x_stride * x_len;
685 } 685 }
686 686
687 int idx_j; 687 octave_idx_type idx_j;
688 688
689 double tmp_max = octave_NaN; 689 double tmp_max = octave_NaN;
690 690
691 for (idx_j = 0; idx_j < x_len; idx_j++) 691 for (idx_j = 0; idx_j < x_len; idx_j++)
692 { 692 {
694 694
695 if (! octave_is_NaN_or_NA (tmp_max)) 695 if (! octave_is_NaN_or_NA (tmp_max))
696 break; 696 break;
697 } 697 }
698 698
699 for (int j = idx_j+1; j < x_len; j++) 699 for (octave_idx_type j = idx_j+1; j < x_len; j++)
700 { 700 {
701 double tmp = elem (j * x_stride + x_offset); 701 double tmp = elem (j * x_stride + x_offset);
702 702
703 if (octave_is_NaN_or_NA (tmp)) 703 if (octave_is_NaN_or_NA (tmp))
704 continue; 704 continue;
717 } 717 }
718 718
719 NDArray 719 NDArray
720 NDArray::min (int dim) const 720 NDArray::min (int dim) const
721 { 721 {
722 ArrayN<int> dummy_idx; 722 ArrayN<octave_idx_type> dummy_idx;
723 return min (dummy_idx, dim); 723 return min (dummy_idx, dim);
724 } 724 }
725 725
726 NDArray 726 NDArray
727 NDArray::min (ArrayN<int>& idx_arg, int dim) const 727 NDArray::min (ArrayN<octave_idx_type>& idx_arg, int dim) const
728 { 728 {
729 dim_vector dv = dims (); 729 dim_vector dv = dims ();
730 dim_vector dr = dims (); 730 dim_vector dr = dims ();
731 731
732 if (dv.numel () == 0 || dim > dv.length () || dim < 0) 732 if (dv.numel () == 0 || dim > dv.length () || dim < 0)
735 dr(dim) = 1; 735 dr(dim) = 1;
736 736
737 NDArray result (dr); 737 NDArray result (dr);
738 idx_arg.resize (dr); 738 idx_arg.resize (dr);
739 739
740 int x_stride = 1; 740 octave_idx_type x_stride = 1;
741 int x_len = dv(dim); 741 octave_idx_type x_len = dv(dim);
742 for (int i = 0; i < dim; i++) 742 for (int i = 0; i < dim; i++)
743 x_stride *= dv(i); 743 x_stride *= dv(i);
744 744
745 for (int i = 0; i < dr.numel (); i++) 745 for (octave_idx_type i = 0; i < dr.numel (); i++)
746 { 746 {
747 int x_offset; 747 octave_idx_type x_offset;
748 if (x_stride == 1) 748 if (x_stride == 1)
749 x_offset = i * x_len; 749 x_offset = i * x_len;
750 else 750 else
751 { 751 {
752 int x_offset2 = 0; 752 octave_idx_type x_offset2 = 0;
753 x_offset = i; 753 x_offset = i;
754 while (x_offset >= x_stride) 754 while (x_offset >= x_stride)
755 { 755 {
756 x_offset -= x_stride; 756 x_offset -= x_stride;
757 x_offset2++; 757 x_offset2++;
758 } 758 }
759 x_offset += x_offset2 * x_stride * x_len; 759 x_offset += x_offset2 * x_stride * x_len;
760 } 760 }
761 761
762 int idx_j; 762 octave_idx_type idx_j;
763 763
764 double tmp_min = octave_NaN; 764 double tmp_min = octave_NaN;
765 765
766 for (idx_j = 0; idx_j < x_len; idx_j++) 766 for (idx_j = 0; idx_j < x_len; idx_j++)
767 { 767 {
769 769
770 if (! octave_is_NaN_or_NA (tmp_min)) 770 if (! octave_is_NaN_or_NA (tmp_min))
771 break; 771 break;
772 } 772 }
773 773
774 for (int j = idx_j+1; j < x_len; j++) 774 for (octave_idx_type j = idx_j+1; j < x_len; j++)
775 { 775 {
776 double tmp = elem (j * x_stride + x_offset); 776 double tmp = elem (j * x_stride + x_offset);
777 777
778 if (octave_is_NaN_or_NA (tmp)) 778 if (octave_is_NaN_or_NA (tmp))
779 continue; 779 continue;
790 790
791 return result; 791 return result;
792 } 792 }
793 793
794 NDArray 794 NDArray
795 NDArray::concat (const NDArray& rb, const Array<int>& ra_idx) 795 NDArray::concat (const NDArray& rb, const Array<octave_idx_type>& ra_idx)
796 { 796 {
797 if (rb.numel () > 0) 797 if (rb.numel () > 0)
798 insert (rb, ra_idx); 798 insert (rb, ra_idx);
799 return *this; 799 return *this;
800 } 800 }
801 801
802 ComplexNDArray 802 ComplexNDArray
803 NDArray::concat (const ComplexNDArray& rb, const Array<int>& ra_idx) 803 NDArray::concat (const ComplexNDArray& rb, const Array<octave_idx_type>& ra_idx)
804 { 804 {
805 ComplexNDArray retval (*this); 805 ComplexNDArray retval (*this);
806 if (rb.numel () > 0) 806 if (rb.numel () > 0)
807 retval.insert (rb, ra_idx); 807 retval.insert (rb, ra_idx);
808 return retval; 808 return retval;
809 } 809 }
810 810
811 charNDArray 811 charNDArray
812 NDArray::concat (const charNDArray& rb, const Array<int>& ra_idx) 812 NDArray::concat (const charNDArray& rb, const Array<octave_idx_type>& ra_idx)
813 { 813 {
814 charNDArray retval (dims ()); 814 charNDArray retval (dims ());
815 int nel = numel (); 815 octave_idx_type nel = numel ();
816 816
817 for (int i = 0; i < nel; i++) 817 for (octave_idx_type i = 0; i < nel; i++)
818 { 818 {
819 double d = elem (i); 819 double d = elem (i);
820 820
821 if (xisnan (d)) 821 if (xisnan (d))
822 { 822 {
824 ("invalid conversion from NaN to character"); 824 ("invalid conversion from NaN to character");
825 return retval; 825 return retval;
826 } 826 }
827 else 827 else
828 { 828 {
829 int ival = NINT (d); 829 octave_idx_type ival = NINTbig (d);
830 830
831 if (ival < 0 || ival > UCHAR_MAX) 831 if (ival < 0 || ival > UCHAR_MAX)
832 // XXX FIXME XXX -- is there something 832 // XXX FIXME XXX -- is there something
833 // better we could do? Should we warn the user? 833 // better we could do? Should we warn the user?
834 ival = 0; 834 ival = 0;
845 } 845 }
846 846
847 NDArray 847 NDArray
848 real (const ComplexNDArray& a) 848 real (const ComplexNDArray& a)
849 { 849 {
850 int a_len = a.length (); 850 octave_idx_type a_len = a.length ();
851 NDArray retval; 851 NDArray retval;
852 if (a_len > 0) 852 if (a_len > 0)
853 retval = NDArray (mx_inline_real_dup (a.data (), a_len), a.dims ()); 853 retval = NDArray (mx_inline_real_dup (a.data (), a_len), a.dims ());
854 return retval; 854 return retval;
855 } 855 }
856 856
857 NDArray 857 NDArray
858 imag (const ComplexNDArray& a) 858 imag (const ComplexNDArray& a)
859 { 859 {
860 int a_len = a.length (); 860 octave_idx_type a_len = a.length ();
861 NDArray retval; 861 NDArray retval;
862 if (a_len > 0) 862 if (a_len > 0)
863 retval = NDArray (mx_inline_imag_dup (a.data (), a_len), a.dims ()); 863 retval = NDArray (mx_inline_imag_dup (a.data (), a_len), a.dims ());
864 return retval; 864 return retval;
865 } 865 }
866 866
867 NDArray& 867 NDArray&
868 NDArray::insert (const NDArray& a, int r, int c) 868 NDArray::insert (const NDArray& a, octave_idx_type r, octave_idx_type c)
869 { 869 {
870 Array<double>::insert (a, r, c); 870 Array<double>::insert (a, r, c);
871 return *this; 871 return *this;
872 } 872 }
873 873
874 NDArray& 874 NDArray&
875 NDArray::insert (const NDArray& a, const Array<int>& ra_idx) 875 NDArray::insert (const NDArray& a, const Array<octave_idx_type>& ra_idx)
876 { 876 {
877 Array<double>::insert (a, ra_idx); 877 Array<double>::insert (a, ra_idx);
878 return *this; 878 return *this;
879 } 879 }
880 880
881 NDArray 881 NDArray
882 NDArray::abs (void) const 882 NDArray::abs (void) const
883 { 883 {
884 NDArray retval (dims ()); 884 NDArray retval (dims ());
885 885
886 int nel = nelem (); 886 octave_idx_type nel = nelem ();
887 887
888 for (int i = 0; i < nel; i++) 888 for (octave_idx_type i = 0; i < nel; i++)
889 retval(i) = fabs (elem (i)); 889 retval(i) = fabs (elem (i));
890 890
891 return retval; 891 return retval;
892 } 892 }
893 893
916 916
917 return retval; 917 return retval;
918 } 918 }
919 919
920 void 920 void
921 NDArray::increment_index (Array<int>& ra_idx, 921 NDArray::increment_index (Array<octave_idx_type>& ra_idx,
922 const dim_vector& dimensions, 922 const dim_vector& dimensions,
923 int start_dimension) 923 int start_dimension)
924 { 924 {
925 ::increment_index (ra_idx, dimensions, start_dimension); 925 ::increment_index (ra_idx, dimensions, start_dimension);
926 } 926 }
927 927
928 int 928 octave_idx_type
929 NDArray::compute_index (Array<int>& ra_idx, 929 NDArray::compute_index (Array<octave_idx_type>& ra_idx,
930 const dim_vector& dimensions) 930 const dim_vector& dimensions)
931 { 931 {
932 return ::compute_index (ra_idx, dimensions); 932 return ::compute_index (ra_idx, dimensions);
933 } 933 }
934 934
935 // This contains no information on the array structure !!! 935 // This contains no information on the array structure !!!
936 std::ostream& 936 std::ostream&
937 operator << (std::ostream& os, const NDArray& a) 937 operator << (std::ostream& os, const NDArray& a)
938 { 938 {
939 int nel = a.nelem (); 939 octave_idx_type nel = a.nelem ();
940 940
941 for (int i = 0; i < nel; i++) 941 for (octave_idx_type i = 0; i < nel; i++)
942 { 942 {
943 os << " "; 943 os << " ";
944 octave_write_double (os, a.elem (i)); 944 octave_write_double (os, a.elem (i));
945 os << "\n"; 945 os << "\n";
946 } 946 }
948 } 948 }
949 949
950 std::istream& 950 std::istream&
951 operator >> (std::istream& is, NDArray& a) 951 operator >> (std::istream& is, NDArray& a)
952 { 952 {
953 int nel = a.nelem (); 953 octave_idx_type nel = a.nelem ();
954 954
955 if (nel < 1 ) 955 if (nel < 1 )
956 is.clear (std::ios::badbit); 956 is.clear (std::ios::badbit);
957 else 957 else
958 { 958 {
959 double tmp; 959 double tmp;
960 for (int i = 0; i < nel; i++) 960 for (octave_idx_type i = 0; i < nel; i++)
961 { 961 {
962 tmp = octave_read_double (is); 962 tmp = octave_read_double (is);
963 if (is) 963 if (is)
964 a.elem (i) = tmp; 964 a.elem (i) = tmp;
965 else 965 else
981 981
982 NDArray 982 NDArray
983 min (double d, const NDArray& m) 983 min (double d, const NDArray& m)
984 { 984 {
985 dim_vector dv = m.dims (); 985 dim_vector dv = m.dims ();
986 int nel = dv.numel (); 986 octave_idx_type nel = dv.numel ();
987 987
988 EMPTY_RETURN_CHECK (NDArray); 988 EMPTY_RETURN_CHECK (NDArray);
989 989
990 NDArray result (dv); 990 NDArray result (dv);
991 991
992 for (int i = 0; i < nel; i++) 992 for (octave_idx_type i = 0; i < nel; i++)
993 { 993 {
994 OCTAVE_QUIT; 994 OCTAVE_QUIT;
995 result (i) = xmin (d, m (i)); 995 result (i) = xmin (d, m (i));
996 } 996 }
997 997
1000 1000
1001 NDArray 1001 NDArray
1002 min (const NDArray& m, double d) 1002 min (const NDArray& m, double d)
1003 { 1003 {
1004 dim_vector dv = m.dims (); 1004 dim_vector dv = m.dims ();
1005 int nel = dv.numel (); 1005 octave_idx_type nel = dv.numel ();
1006 1006
1007 EMPTY_RETURN_CHECK (NDArray); 1007 EMPTY_RETURN_CHECK (NDArray);
1008 1008
1009 NDArray result (dv); 1009 NDArray result (dv);
1010 1010
1011 for (int i = 0; i < nel; i++) 1011 for (octave_idx_type i = 0; i < nel; i++)
1012 { 1012 {
1013 OCTAVE_QUIT; 1013 OCTAVE_QUIT;
1014 result (i) = xmin (d, m (i)); 1014 result (i) = xmin (d, m (i));
1015 } 1015 }
1016 1016
1019 1019
1020 NDArray 1020 NDArray
1021 min (const NDArray& a, const NDArray& b) 1021 min (const NDArray& a, const NDArray& b)
1022 { 1022 {
1023 dim_vector dv = a.dims (); 1023 dim_vector dv = a.dims ();
1024 int nel = dv.numel (); 1024 octave_idx_type nel = dv.numel ();
1025 1025
1026 if (dv != b.dims ()) 1026 if (dv != b.dims ())
1027 { 1027 {
1028 (*current_liboctave_error_handler) 1028 (*current_liboctave_error_handler)
1029 ("two-arg min expecting args of same size"); 1029 ("two-arg min expecting args of same size");
1032 1032
1033 EMPTY_RETURN_CHECK (NDArray); 1033 EMPTY_RETURN_CHECK (NDArray);
1034 1034
1035 NDArray result (dv); 1035 NDArray result (dv);
1036 1036
1037 for (int i = 0; i < nel; i++) 1037 for (octave_idx_type i = 0; i < nel; i++)
1038 { 1038 {
1039 OCTAVE_QUIT; 1039 OCTAVE_QUIT;
1040 result (i) = xmin (a (i), b (i)); 1040 result (i) = xmin (a (i), b (i));
1041 } 1041 }
1042 1042
1045 1045
1046 NDArray 1046 NDArray
1047 max (double d, const NDArray& m) 1047 max (double d, const NDArray& m)
1048 { 1048 {
1049 dim_vector dv = m.dims (); 1049 dim_vector dv = m.dims ();
1050 int nel = dv.numel (); 1050 octave_idx_type nel = dv.numel ();
1051 1051
1052 EMPTY_RETURN_CHECK (NDArray); 1052 EMPTY_RETURN_CHECK (NDArray);
1053 1053
1054 NDArray result (dv); 1054 NDArray result (dv);
1055 1055
1056 for (int i = 0; i < nel; i++) 1056 for (octave_idx_type i = 0; i < nel; i++)
1057 { 1057 {
1058 OCTAVE_QUIT; 1058 OCTAVE_QUIT;
1059 result (i) = xmax (d, m (i)); 1059 result (i) = xmax (d, m (i));
1060 } 1060 }
1061 1061
1064 1064
1065 NDArray 1065 NDArray
1066 max (const NDArray& m, double d) 1066 max (const NDArray& m, double d)
1067 { 1067 {
1068 dim_vector dv = m.dims (); 1068 dim_vector dv = m.dims ();
1069 int nel = dv.numel (); 1069 octave_idx_type nel = dv.numel ();
1070 1070
1071 EMPTY_RETURN_CHECK (NDArray); 1071 EMPTY_RETURN_CHECK (NDArray);
1072 1072
1073 NDArray result (dv); 1073 NDArray result (dv);
1074 1074
1075 for (int i = 0; i < nel; i++) 1075 for (octave_idx_type i = 0; i < nel; i++)
1076 { 1076 {
1077 OCTAVE_QUIT; 1077 OCTAVE_QUIT;
1078 result (i) = xmax (d, m (i)); 1078 result (i) = xmax (d, m (i));
1079 } 1079 }
1080 1080
1083 1083
1084 NDArray 1084 NDArray
1085 max (const NDArray& a, const NDArray& b) 1085 max (const NDArray& a, const NDArray& b)
1086 { 1086 {
1087 dim_vector dv = a.dims (); 1087 dim_vector dv = a.dims ();
1088 int nel = dv.numel (); 1088 octave_idx_type nel = dv.numel ();
1089 1089
1090 if (dv != b.dims ()) 1090 if (dv != b.dims ())
1091 { 1091 {
1092 (*current_liboctave_error_handler) 1092 (*current_liboctave_error_handler)
1093 ("two-arg max expecting args of same size"); 1093 ("two-arg max expecting args of same size");
1096 1096
1097 EMPTY_RETURN_CHECK (NDArray); 1097 EMPTY_RETURN_CHECK (NDArray);
1098 1098
1099 NDArray result (dv); 1099 NDArray result (dv);
1100 1100
1101 for (int i = 0; i < nel; i++) 1101 for (octave_idx_type i = 0; i < nel; i++)
1102 { 1102 {
1103 OCTAVE_QUIT; 1103 OCTAVE_QUIT;
1104 result (i) = xmax (a (i), b (i)); 1104 result (i) = xmax (a (i), b (i));
1105 } 1105 }
1106 1106