Mercurial > hg > octave-nkf

comparison liboctave/fCNDArray.cc @ 10314:07ebe522dac2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
untabify liboctave C++ sources
author John W. Eaton <jwe@octave.org>
date Thu, 11 Feb 2010 12:23:32 -0500
parents 4c0cdbe0acca
children 83fa590b8a09
comparison
equal deleted inserted replaced
10313:f3b65e1ae355 10314:07ebe522dac2
79 FloatComplex *out (retval.fortran_vec ()); 79 FloatComplex *out (retval.fortran_vec ());
80 80
81 // Need to be careful here about the distance between fft's 81 // Need to be careful here about the distance between fft's
82 for (octave_idx_type k = 0; k < nloop; k++) 82 for (octave_idx_type k = 0; k < nloop; k++)
83 octave_fftw::fft (in + k * stride * n, out + k * stride * n, 83 octave_fftw::fft (in + k * stride * n, out + k * stride * n,
84 n, howmany, stride, dist); 84 n, howmany, stride, dist);
85 85
86 return retval; 86 return retval;
87 } 87 }
88 88
89 FloatComplexNDArray 89 FloatComplexNDArray
110 FloatComplex *out (retval.fortran_vec ()); 110 FloatComplex *out (retval.fortran_vec ());
111 111
112 // Need to be careful here about the distance between fft's 112 // Need to be careful here about the distance between fft's
113 for (octave_idx_type k = 0; k < nloop; k++) 113 for (octave_idx_type k = 0; k < nloop; k++)
114 octave_fftw::ifft (in + k * stride * n, out + k * stride * n, 114 octave_fftw::ifft (in + k * stride * n, out + k * stride * n,
115 n, howmany, stride, dist); 115 n, howmany, stride, dist);
116 116
117 return retval; 117 return retval;
118 } 118 }
119 119
120 FloatComplexNDArray 120 FloatComplexNDArray
230 F77_FUNC (cffti, CFFTI) (npts, pwsave); 230 F77_FUNC (cffti, CFFTI) (npts, pwsave);
231 231
232 for (octave_idx_type k = 0; k < nloop; k++) 232 for (octave_idx_type k = 0; k < nloop; k++)
233 { 233 {
234 for (octave_idx_type j = 0; j < howmany; j++) 234 for (octave_idx_type j = 0; j < howmany; j++)
235 { 235 {
236 octave_quit (); 236 octave_quit ();
237 237
238 for (octave_idx_type i = 0; i < npts; i++) 238 for (octave_idx_type i = 0; i < npts; i++)
239 tmp[i] = elem((i + k*npts)*stride + j*dist); 239 tmp[i] = elem((i + k*npts)*stride + j*dist);
240 240
241 F77_FUNC (cfftf, CFFTF) (npts, tmp, pwsave); 241 F77_FUNC (cfftf, CFFTF) (npts, tmp, pwsave);
242 242
243 for (octave_idx_type i = 0; i < npts; i++) 243 for (octave_idx_type i = 0; i < npts; i++)
244 retval ((i + k*npts)*stride + j*dist) = tmp[i]; 244 retval ((i + k*npts)*stride + j*dist) = tmp[i];
245 } 245 }
246 } 246 }
247 247
248 return retval; 248 return retval;
249 } 249 }
250 250
277 F77_FUNC (cffti, CFFTI) (npts, pwsave); 277 F77_FUNC (cffti, CFFTI) (npts, pwsave);
278 278
279 for (octave_idx_type k = 0; k < nloop; k++) 279 for (octave_idx_type k = 0; k < nloop; k++)
280 { 280 {
281 for (octave_idx_type j = 0; j < howmany; j++) 281 for (octave_idx_type j = 0; j < howmany; j++)
282 { 282 {
283 octave_quit (); 283 octave_quit ();
284 284
285 for (octave_idx_type i = 0; i < npts; i++) 285 for (octave_idx_type i = 0; i < npts; i++)
286 tmp[i] = elem((i + k*npts)*stride + j*dist); 286 tmp[i] = elem((i + k*npts)*stride + j*dist);
287 287
288 F77_FUNC (cfftb, CFFTB) (npts, tmp, pwsave); 288 F77_FUNC (cfftb, CFFTB) (npts, tmp, pwsave);
289 289
290 for (octave_idx_type i = 0; i < npts; i++) 290 for (octave_idx_type i = 0; i < npts; i++)
291 retval ((i + k*npts)*stride + j*dist) = tmp[i] / 291 retval ((i + k*npts)*stride + j*dist) = tmp[i] /
292 static_cast<float> (npts); 292 static_cast<float> (npts);
293 } 293 }
294 } 294 }
295 295
296 return retval; 296 return retval;
297 } 297 }
298 298
314 Array<FloatComplex> row (npts); 314 Array<FloatComplex> row (npts);
315 FloatComplex *prow = row.fortran_vec (); 315 FloatComplex *prow = row.fortran_vec ();
316 316
317 octave_idx_type howmany = numel () / npts; 317 octave_idx_type howmany = numel () / npts;
318 howmany = (stride == 1 ? howmany : 318 howmany = (stride == 1 ? howmany :
319 (howmany > stride ? stride : howmany)); 319 (howmany > stride ? stride : howmany));
320 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride); 320 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride);
321 octave_idx_type dist = (stride == 1 ? npts : 1); 321 octave_idx_type dist = (stride == 1 ? npts : 1);
322 322
323 F77_FUNC (cffti, CFFTI) (npts, pwsave); 323 F77_FUNC (cffti, CFFTI) (npts, pwsave);
324 324
325 for (octave_idx_type k = 0; k < nloop; k++) 325 for (octave_idx_type k = 0; k < nloop; k++)
326 { 326 {
327 for (octave_idx_type j = 0; j < howmany; j++) 327 for (octave_idx_type j = 0; j < howmany; j++)
328 { 328 {
329 octave_quit (); 329 octave_quit ();
330 330
331 for (octave_idx_type l = 0; l < npts; l++) 331 for (octave_idx_type l = 0; l < npts; l++)
332 prow[l] = retval ((l + k*npts)*stride + j*dist); 332 prow[l] = retval ((l + k*npts)*stride + j*dist);
333 333
334 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave); 334 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
335 335
336 for (octave_idx_type l = 0; l < npts; l++) 336 for (octave_idx_type l = 0; l < npts; l++)
337 retval ((l + k*npts)*stride + j*dist) = prow[l]; 337 retval ((l + k*npts)*stride + j*dist) = prow[l];
338 } 338 }
339 } 339 }
340 340
341 stride *= dv2(i); 341 stride *= dv2(i);
342 } 342 }
343 343
344 return retval; 344 return retval;
362 Array<FloatComplex> row (npts); 362 Array<FloatComplex> row (npts);
363 FloatComplex *prow = row.fortran_vec (); 363 FloatComplex *prow = row.fortran_vec ();
364 364
365 octave_idx_type howmany = numel () / npts; 365 octave_idx_type howmany = numel () / npts;
366 howmany = (stride == 1 ? howmany : 366 howmany = (stride == 1 ? howmany :
367 (howmany > stride ? stride : howmany)); 367 (howmany > stride ? stride : howmany));
368 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride); 368 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride);
369 octave_idx_type dist = (stride == 1 ? npts : 1); 369 octave_idx_type dist = (stride == 1 ? npts : 1);
370 370
371 F77_FUNC (cffti, CFFTI) (npts, pwsave); 371 F77_FUNC (cffti, CFFTI) (npts, pwsave);
372 372
373 for (octave_idx_type k = 0; k < nloop; k++) 373 for (octave_idx_type k = 0; k < nloop; k++)
374 { 374 {
375 for (octave_idx_type j = 0; j < howmany; j++) 375 for (octave_idx_type j = 0; j < howmany; j++)
376 { 376 {
377 octave_quit (); 377 octave_quit ();
378 378
379 for (octave_idx_type l = 0; l < npts; l++) 379 for (octave_idx_type l = 0; l < npts; l++)
380 prow[l] = retval ((l + k*npts)*stride + j*dist); 380 prow[l] = retval ((l + k*npts)*stride + j*dist);
381 381
382 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave); 382 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
383 383
384 for (octave_idx_type l = 0; l < npts; l++) 384 for (octave_idx_type l = 0; l < npts; l++)
385 retval ((l + k*npts)*stride + j*dist) = prow[l] / 385 retval ((l + k*npts)*stride + j*dist) = prow[l] /
386 static_cast<float> (npts); 386 static_cast<float> (npts);
387 } 387 }
388 } 388 }
389 389
390 stride *= dv2(i); 390 stride *= dv2(i);
391 } 391 }
392 392
393 return retval; 393 return retval;
410 Array<FloatComplex> row (npts); 410 Array<FloatComplex> row (npts);
411 FloatComplex *prow = row.fortran_vec (); 411 FloatComplex *prow = row.fortran_vec ();
412 412
413 octave_idx_type howmany = numel () / npts; 413 octave_idx_type howmany = numel () / npts;
414 howmany = (stride == 1 ? howmany : 414 howmany = (stride == 1 ? howmany :
415 (howmany > stride ? stride : howmany)); 415 (howmany > stride ? stride : howmany));
416 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride); 416 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride);
417 octave_idx_type dist = (stride == 1 ? npts : 1); 417 octave_idx_type dist = (stride == 1 ? npts : 1);
418 418
419 F77_FUNC (cffti, CFFTI) (npts, pwsave); 419 F77_FUNC (cffti, CFFTI) (npts, pwsave);
420 420
421 for (octave_idx_type k = 0; k < nloop; k++) 421 for (octave_idx_type k = 0; k < nloop; k++)
422 { 422 {
423 for (octave_idx_type j = 0; j < howmany; j++) 423 for (octave_idx_type j = 0; j < howmany; j++)
424 { 424 {
425 octave_quit (); 425 octave_quit ();
426 426
427 for (octave_idx_type l = 0; l < npts; l++) 427 for (octave_idx_type l = 0; l < npts; l++)
428 prow[l] = retval ((l + k*npts)*stride + j*dist); 428 prow[l] = retval ((l + k*npts)*stride + j*dist);
429 429
430 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave); 430 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
431 431
432 for (octave_idx_type l = 0; l < npts; l++) 432 for (octave_idx_type l = 0; l < npts; l++)
433 retval ((l + k*npts)*stride + j*dist) = prow[l]; 433 retval ((l + k*npts)*stride + j*dist) = prow[l];
434 } 434 }
435 } 435 }
436 436
437 stride *= dv(i); 437 stride *= dv(i);
438 } 438 }
439 439
440 return retval; 440 return retval;
457 Array<FloatComplex> row (npts); 457 Array<FloatComplex> row (npts);
458 FloatComplex *prow = row.fortran_vec (); 458 FloatComplex *prow = row.fortran_vec ();
459 459
460 octave_idx_type howmany = numel () / npts; 460 octave_idx_type howmany = numel () / npts;
461 howmany = (stride == 1 ? howmany : 461 howmany = (stride == 1 ? howmany :
462 (howmany > stride ? stride : howmany)); 462 (howmany > stride ? stride : howmany));
463 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride); 463 octave_idx_type nloop = (stride == 1 ? 1 : numel () / npts / stride);
464 octave_idx_type dist = (stride == 1 ? npts : 1); 464 octave_idx_type dist = (stride == 1 ? npts : 1);
465 465
466 F77_FUNC (cffti, CFFTI) (npts, pwsave); 466 F77_FUNC (cffti, CFFTI) (npts, pwsave);
467 467
468 for (octave_idx_type k = 0; k < nloop; k++) 468 for (octave_idx_type k = 0; k < nloop; k++)
469 { 469 {
470 for (octave_idx_type j = 0; j < howmany; j++) 470 for (octave_idx_type j = 0; j < howmany; j++)
471 { 471 {
472 octave_quit (); 472 octave_quit ();
473 473
474 for (octave_idx_type l = 0; l < npts; l++) 474 for (octave_idx_type l = 0; l < npts; l++)
475 prow[l] = retval ((l + k*npts)*stride + j*dist); 475 prow[l] = retval ((l + k*npts)*stride + j*dist);
476 476
477 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave); 477 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
478 478
479 for (octave_idx_type l = 0; l < npts; l++) 479 for (octave_idx_type l = 0; l < npts; l++)
480 retval ((l + k*npts)*stride + j*dist) = prow[l] / 480 retval ((l + k*npts)*stride + j*dist) = prow[l] /
481 static_cast<float> (npts); 481 static_cast<float> (npts);
482 } 482 }
483 } 483 }
484 484
485 stride *= dv(i); 485 stride *= dv(i);
486 } 486 }
487 487
488 return retval; 488 return retval;
507 507
508 for (octave_idx_type i = 0; i < nel; i++) 508 for (octave_idx_type i = 0; i < nel; i++)
509 { 509 {
510 FloatComplex val = elem (i); 510 FloatComplex val = elem (i);
511 if (xisnan (val)) 511 if (xisnan (val))
512 return true; 512 return true;
513 } 513 }
514 return false; 514 return false;
515 } 515 }
516 516
517 bool 517 bool
521 521
522 for (octave_idx_type i = 0; i < nel; i++) 522 for (octave_idx_type i = 0; i < nel; i++)
523 { 523 {
524 FloatComplex val = elem (i); 524 FloatComplex val = elem (i);
525 if (xisinf (val) || xisnan (val)) 525 if (xisinf (val) || xisnan (val))
526 return true; 526 return true;
527 } 527 }
528 return false; 528 return false;
529 } 529 }
530 530
531 // Return true if no elements have imaginary components. 531 // Return true if no elements have imaginary components.
554 554
555 max_val = r_val; 555 max_val = r_val;
556 min_val = r_val; 556 min_val = r_val;
557 557
558 if (i_val > max_val) 558 if (i_val > max_val)
559 max_val = i_val; 559 max_val = i_val;
560 560
561 if (i_val < max_val) 561 if (i_val < max_val)
562 min_val = i_val; 562 min_val = i_val;
563 } 563 }
564 else 564 else
565 return false; 565 return false;
566 566
567 for (octave_idx_type i = 0; i < nel; i++) 567 for (octave_idx_type i = 0; i < nel; i++)
570 570
571 float r_val = std::real (val); 571 float r_val = std::real (val);
572 float i_val = std::imag (val); 572 float i_val = std::imag (val);
573 573
574 if (r_val > max_val) 574 if (r_val > max_val)
575 max_val = r_val; 575 max_val = r_val;
576 576
577 if (i_val > max_val) 577 if (i_val > max_val)
578 max_val = i_val; 578 max_val = i_val;
579 579
580 if (r_val < min_val) 580 if (r_val < min_val)
581 min_val = r_val; 581 min_val = r_val;
582 582
583 if (i_val < min_val) 583 if (i_val < min_val)
584 min_val = i_val; 584 min_val = i_val;
585 585
586 if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val) 586 if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val)
587 return false; 587 return false;
588 } 588 }
589 589
590 return true; 590 return true;
591 } 591 }
592 592
601 601
602 float r_val = std::real (val); 602 float r_val = std::real (val);
603 float i_val = std::imag (val); 603 float i_val = std::imag (val);
604 604
605 if ((! (xisnan (r_val) || xisinf (r_val)) 605 if ((! (xisnan (r_val) || xisinf (r_val))
606 && fabs (r_val) > FLT_MAX) 606 && fabs (r_val) > FLT_MAX)
607 || (! (xisnan (i_val) || xisinf (i_val)) 607 || (! (xisnan (i_val) || xisinf (i_val))
608 && fabs (i_val) > FLT_MAX)) 608 && fabs (i_val) > FLT_MAX))
609 return true; 609 return true;
610 } 610 }
611 611
612 return false; 612 return false;
613 } 613 }
614 614
785 785
786 a_ra_idx.elem (0) = r; 786 a_ra_idx.elem (0) = r;
787 a_ra_idx.elem (1) = c; 787 a_ra_idx.elem (1) = c;
788 788
789 for (int i = 0; i < n; i++) 789 for (int i = 0; i < n; i++)
790 { 790 {
791 if (a_ra_idx (i) < 0 || (a_ra_idx (i) + a_dv (i)) > dimensions (i)) 791 if (a_ra_idx (i) < 0 || (a_ra_idx (i) + a_dv (i)) > dimensions (i))
792 { 792 {
793 (*current_liboctave_error_handler) 793 (*current_liboctave_error_handler)
794 ("Array<T>::insert: range error for insert"); 794 ("Array<T>::insert: range error for insert");
795 return *this; 795 return *this;
796 } 796 }
797 } 797 }
798 798
799 a_ra_idx.elem (0) = 0; 799 a_ra_idx.elem (0) = 0;
800 a_ra_idx.elem (1) = 0; 800 a_ra_idx.elem (1) = 0;
801 801
802 octave_idx_type n_elt = a.numel (); 802 octave_idx_type n_elt = a.numel ();
803 803
804 // IS make_unique () NECCESSARY HERE?? 804 // IS make_unique () NECCESSARY HERE??
805 805
806 for (octave_idx_type i = 0; i < n_elt; i++) 806 for (octave_idx_type i = 0; i < n_elt; i++)
807 { 807 {
808 Array<octave_idx_type> ra_idx = a_ra_idx; 808 Array<octave_idx_type> ra_idx = a_ra_idx;
809 809
810 ra_idx.elem (0) = a_ra_idx (0) + r; 810 ra_idx.elem (0) = a_ra_idx (0) + r;
811 ra_idx.elem (1) = a_ra_idx (1) + c; 811 ra_idx.elem (1) = a_ra_idx (1) + c;
812 812
813 elem (ra_idx) = a.elem (a_ra_idx); 813 elem (ra_idx) = a.elem (a_ra_idx);
814 814
815 increment_index (a_ra_idx, a_dv); 815 increment_index (a_ra_idx, a_dv);
816 } 816 }
817 } 817 }
818 else 818 else
819 (*current_liboctave_error_handler) 819 (*current_liboctave_error_handler)
820 ("Array<T>::insert: invalid indexing operation"); 820 ("Array<T>::insert: invalid indexing operation");
821 821
850 return retval; 850 return retval;
851 } 851 }
852 852
853 void 853 void
854 FloatComplexNDArray::increment_index (Array<octave_idx_type>& ra_idx, 854 FloatComplexNDArray::increment_index (Array<octave_idx_type>& ra_idx,
855 const dim_vector& dimensions, 855 const dim_vector& dimensions,
856 int start_dimension) 856 int start_dimension)
857 { 857 {
858 ::increment_index (ra_idx, dimensions, start_dimension); 858 ::increment_index (ra_idx, dimensions, start_dimension);
859 } 859 }
860 860
861 octave_idx_type 861 octave_idx_type
862 FloatComplexNDArray::compute_index (Array<octave_idx_type>& ra_idx, 862 FloatComplexNDArray::compute_index (Array<octave_idx_type>& ra_idx,
863 const dim_vector& dimensions) 863 const dim_vector& dimensions)
864 { 864 {
865 return ::compute_index (ra_idx, dimensions); 865 return ::compute_index (ra_idx, dimensions);
866 } 866 }
867 867
868 FloatComplexNDArray 868 FloatComplexNDArray
893 893
894 if (nel > 0) 894 if (nel > 0)
895 { 895 {
896 FloatComplex tmp; 896 FloatComplex tmp;
897 for (octave_idx_type i = 0; i < nel; i++) 897 for (octave_idx_type i = 0; i < nel; i++)
898 { 898 {
899 tmp = octave_read_value<FloatComplex> (is); 899 tmp = octave_read_value<FloatComplex> (is);
900 if (is) 900 if (is)
901 a.elem (i) = tmp; 901 a.elem (i) = tmp;
902 else 902 else
903 goto done; 903 goto done;
904 } 904 }
905 } 905 }
906 906
907 done: 907 done:
908 908
909 return is; 909 return is;