Mercurial > forge

comparison extra/NaN/src/sumskipnan_mex.cpp @ 12691:6d6285a2a633 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
author schloegl
date Sat, 12 Sep 2015 14:16:39 +0000
parents f26b1170ea90
children 13815b367946
comparison
equal deleted inserted replaced
12690:16bc2657b1f1 12691:6d6285a2a633
49 49
50 50
51 #include <math.h> 51 #include <math.h>
52 #include <stdint.h> 52 #include <stdint.h>
53 #include "mex.h" 53 #include "mex.h"
54
55 /*
56 math.h has isnan() defined for all sizes of floating point numbers,
57 but c++ assumes isnan(double), causing possible conversions for float and long double
58 */
59 #define ISNAN(a) (a!=a)
60
54 61
55 inline void __sumskipnan2w__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W); 62 inline void __sumskipnan2w__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W);
56 inline void __sumskipnan3w__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W); 63 inline void __sumskipnan3w__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W);
57 inline void __sumskipnan2wr__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W); 64 inline void __sumskipnan2wr__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W);
58 inline void __sumskipnan3wr__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W); 65 inline void __sumskipnan3wr__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W);
299 // minimize cache misses 306 // minimize cache misses
300 ix2 = ix0; // index for output 307 ix2 = ix0; // index for output
301 // Inner LOOP: along dimensions < DIM 308 // Inner LOOP: along dimensions < DIM
302 if (W) do { 309 if (W) do {
303 long double x = *LInput; 310 long double x = *LInput;
304 if (!isnan(x)) { 311 if (!ISNAN(x)) {
305 LongOutputSum[ix2] += W[j]*x; 312 LongOutputSum[ix2] += W[j]*x;
306 } 313 }
307 #ifndef NO_FLAG 314 #ifndef NO_FLAG
308 else 315 else
309 flag_isNaN = 1; 316 flag_isNaN = 1;
311 LInput++; 318 LInput++;
312 ix2++; 319 ix2++;
313 } while (ix2 != (l+1)*D1); 320 } while (ix2 != (l+1)*D1);
314 else do { 321 else do {
315 long double x = *LInput; 322 long double x = *LInput;
316 if (!isnan(x)) { 323 if (!ISNAN(x)) {
317 LongOutputSum[ix2] += x; 324 LongOutputSum[ix2] += x;
318 } 325 }
319 #ifndef NO_FLAG 326 #ifndef NO_FLAG
320 else 327 else
321 flag_isNaN = 1; 328 flag_isNaN = 1;
341 // minimize cache misses 348 // minimize cache misses
342 ix2 = ix0; // index for output 349 ix2 = ix0; // index for output
343 // Inner LOOP: along dimensions < DIM 350 // Inner LOOP: along dimensions < DIM
344 if (W) do { 351 if (W) do {
345 long double x = *LInput; 352 long double x = *LInput;
346 if (!isnan(x)) { 353 if (!ISNAN(x)) {
347 LongOutputCount[ix2] += W[j]; 354 LongOutputCount[ix2] += W[j];
348 LongOutputSum[ix2] += W[j]*x; 355 LongOutputSum[ix2] += W[j]*x;
349 } 356 }
350 #ifndef NO_FLAG 357 #ifndef NO_FLAG
351 else 358 else
354 LInput++; 361 LInput++;
355 ix2++; 362 ix2++;
356 } while (ix2 != (l+1)*D1); 363 } while (ix2 != (l+1)*D1);
357 else do { 364 else do {
358 long double x = *LInput; 365 long double x = *LInput;
359 if (!isnan(x)) { 366 if (!ISNAN(x)) {
360 LongOutputCount[ix2] += 1.0; 367 LongOutputCount[ix2] += 1.0;
361 LongOutputSum[ix2] += x; 368 LongOutputSum[ix2] += x;
362 } 369 }
363 #ifndef NO_FLAG 370 #ifndef NO_FLAG
364 else 371 else
386 // minimize cache misses 393 // minimize cache misses
387 ix2 = ix0; // index for output 394 ix2 = ix0; // index for output
388 // Inner LOOP: along dimensions < DIM 395 // Inner LOOP: along dimensions < DIM
389 if (W) do { 396 if (W) do {
390 long double x = *LInput; 397 long double x = *LInput;
391 if (!isnan(x)) { 398 if (!ISNAN(x)) {
392 LongOutputCount[ix2] += W[j]; 399 LongOutputCount[ix2] += W[j];
393 long double t = W[j]*x; 400 long double t = W[j]*x;
394 LongOutputSum[ix2] += t; 401 LongOutputSum[ix2] += t;
395 LongOutputSum2[ix2] += x*t; 402 LongOutputSum2[ix2] += x*t;
396 } 403 }
401 LInput++; 408 LInput++;
402 ix2++; 409 ix2++;
403 } while (ix2 != (l+1)*D1); 410 } while (ix2 != (l+1)*D1);
404 else do { 411 else do {
405 long double x = *LInput; 412 long double x = *LInput;
406 if (!isnan(x)) { 413 if (!ISNAN(x)) {
407 LongOutputCount[ix2] += 1.0; 414 LongOutputCount[ix2] += 1.0;
408 LongOutputSum[ix2] += x; 415 LongOutputSum[ix2] += x;
409 LongOutputSum2[ix2] += x*x; 416 LongOutputSum2[ix2] += x*x;
410 } 417 }
411 #ifndef NO_FLAG 418 #ifndef NO_FLAG
484 if (W) { 491 if (W) {
485 // with weight vector 492 // with weight vector
486 long double count = 0.0; 493 long double count = 0.0;
487 do { 494 do {
488 long double x = *data; 495 long double x = *data;
489 if (!isnan(x)) 496 if (!ISNAN(x))
490 { 497 {
491 count += *W; 498 count += *W;
492 sum += *W*x; 499 sum += *W*x;
493 } 500 }
494 #ifndef NO_FLAG 501 #ifndef NO_FLAG
504 } else { 511 } else {
505 // w/o weight vector 512 // w/o weight vector
506 size_t countI = 0; 513 size_t countI = 0;
507 do { 514 do {
508 long double x = *data; 515 long double x = *data;
509 if (!isnan(x)) 516 if (!ISNAN(x))
510 { 517 {
511 countI++; 518 countI++;
512 sum += x; 519 sum += x;
513 } 520 }
514 #ifndef NO_FLAG 521 #ifndef NO_FLAG
540 if (W) { 547 if (W) {
541 // with weight vector 548 // with weight vector
542 long double count = 0.0; 549 long double count = 0.0;
543 do { 550 do {
544 long double x = *data; 551 long double x = *data;
545 if (!isnan(x)) { 552 if (!ISNAN(x)) {
546 count += *W; 553 count += *W;
547 long double t = *W*x; 554 long double t = *W*x;
548 sum += t; 555 sum += t;
549 msq += x*t; 556 msq += x*t;
550 } 557 }
560 } else { 567 } else {
561 // w/o weight vector 568 // w/o weight vector
562 size_t countI = 0; 569 size_t countI = 0;
563 do { 570 do {
564 long double x = *data; 571 long double x = *data;
565 if (!isnan(x)) { 572 if (!ISNAN(x)) {
566 countI++; 573 countI++;
567 sum += x; 574 sum += x;
568 msq += x*x; 575 msq += x*x;
569 } 576 }
570 #ifndef NO_FLAG 577 #ifndef NO_FLAG
594 if (W) { 601 if (W) {
595 // with weight vector 602 // with weight vector
596 double count = 0.0; 603 double count = 0.0;
597 do { 604 do {
598 double x = *data; 605 double x = *data;
599 if (!isnan(x)) 606 if (!ISNAN(x))
600 { 607 {
601 count += *W; 608 count += *W;
602 sum += *W*x; 609 sum += *W*x;
603 } 610 }
604 #ifndef NO_FLAG 611 #ifndef NO_FLAG
614 } else { 621 } else {
615 // w/o weight vector 622 // w/o weight vector
616 size_t countI = 0; 623 size_t countI = 0;
617 do { 624 do {
618 double x = *data; 625 double x = *data;
619 if (!isnan(x)) 626 if (!ISNAN(x))
620 { 627 {
621 countI++; 628 countI++;
622 sum += x; 629 sum += x;
623 } 630 }
624 #ifndef NO_FLAG 631 #ifndef NO_FLAG
650 if (W) { 657 if (W) {
651 // with weight vector 658 // with weight vector
652 double count = 0.0; 659 double count = 0.0;
653 do { 660 do {
654 double x = *data; 661 double x = *data;
655 if (!isnan(x)) { 662 if (!ISNAN(x)) {
656 count += *W; 663 count += *W;
657 double t = *W*x; 664 double t = *W*x;
658 sum += t; 665 sum += t;
659 msq += x*t; 666 msq += x*t;
660 } 667 }
670 } else { 677 } else {
671 // w/o weight vector 678 // w/o weight vector
672 size_t countI = 0; 679 size_t countI = 0;
673 do { 680 do {
674 double x = *data; 681 double x = *data;
675 if (!isnan(x)) { 682 if (!ISNAN(x)) {
676 countI++; 683 countI++;
677 sum += x; 684 sum += x;
678 msq += x*x; 685 msq += x*x;
679 } 686 }
680 #ifndef NO_FLAG 687 #ifndef NO_FLAG
718 long double count = 0.0; 725 long double count = 0.0;
719 long double rc=0.0, rn=0.0; 726 long double rc=0.0, rn=0.0;
720 do { 727 do {
721 long double x = *data; 728 long double x = *data;
722 long double t,y; 729 long double t,y;
723 if (!isnan(x)) 730 if (!ISNAN(x))
724 { 731 {
725 //count += *W; [1] 732 //count += *W; [1]
726 y = *W-rn; 733 y = *W-rn;
727 t = count+y; 734 t = count+y;
728 rn= (t-count)-y; 735 rn= (t-count)-y;
749 size_t countI = 0; 756 size_t countI = 0;
750 long double rc=0.0; 757 long double rc=0.0;
751 do { 758 do {
752 long double x = *data; 759 long double x = *data;
753 long double t,y; 760 long double t,y;
754 if (!isnan(x)) 761 if (!ISNAN(x))
755 { 762 {
756 countI++; 763 countI++;
757 // sum += x; [1] 764 // sum += x; [1]
758 y = x-rc; 765 y = x-rc;
759 t = sum+y; 766 t = sum+y;
791 long double count = 0.0; 798 long double count = 0.0;
792 long double rc=0.0, rn=0.0, rq=0.0; 799 long double rc=0.0, rn=0.0, rq=0.0;
793 do { 800 do {
794 long double x = *data; 801 long double x = *data;
795 long double t,y; 802 long double t,y;
796 if (!isnan(x)) { 803 if (!ISNAN(x)) {
797 //count += *W; [1] 804 //count += *W; [1]
798 y = *W-rn; 805 y = *W-rn;
799 t = count+y; 806 t = count+y;
800 rn= (t-count)-y; 807 rn= (t-count)-y;
801 count= t; 808 count= t;
827 size_t countI = 0; 834 size_t countI = 0;
828 long double rc=0.0, rq=0.0; 835 long double rc=0.0, rq=0.0;
829 do { 836 do {
830 long double x = *data; 837 long double x = *data;
831 long double t,y; 838 long double t,y;
832 if (!isnan(x)) { 839 if (!ISNAN(x)) {
833 countI++; 840 countI++;
834 //sum += x; [1] 841 //sum += x; [1]
835 y = x-rc; 842 y = x-rc;
836 t = sum+y; 843 t = sum+y;
837 rc= (t-sum)-y; 844 rc= (t-sum)-y;
872 double count = 0.0; 879 double count = 0.0;
873 double rc=0.0, rn=0.0; 880 double rc=0.0, rn=0.0;
874 do { 881 do {
875 double x = *data; 882 double x = *data;
876 double t,y; 883 double t,y;
877 if (!isnan(x)) 884 if (!ISNAN(x))
878 { 885 {
879 //count += *W; [1] 886 //count += *W; [1]
880 y = *W-rn; 887 y = *W-rn;
881 t = count+y; 888 t = count+y;
882 rn= (t-count)-y; 889 rn= (t-count)-y;
903 size_t countI = 0; 910 size_t countI = 0;
904 double rc=0.0; 911 double rc=0.0;
905 do { 912 do {
906 double x = *data; 913 double x = *data;
907 double t,y; 914 double t,y;
908 if (!isnan(x)) 915 if (!ISNAN(x))
909 { 916 {
910 countI++; 917 countI++;
911 // sum += x; [1] 918 // sum += x; [1]
912 y = x-rc; 919 y = x-rc;
913 t = sum+y; 920 t = sum+y;
945 double count = 0.0; 952 double count = 0.0;
946 double rc=0.0, rn=0.0, rq=0.0; 953 double rc=0.0, rn=0.0, rq=0.0;
947 do { 954 do {
948 double x = *data; 955 double x = *data;
949 double t,y; 956 double t,y;
950 if (!isnan(x)) { 957 if (!ISNAN(x)) {
951 //count += *W; [1] 958 //count += *W; [1]
952 y = *W-rn; 959 y = *W-rn;
953 t = count+y; 960 t = count+y;
954 rn= (t-count)-y; 961 rn= (t-count)-y;
955 count= t; 962 count= t;
981 size_t countI = 0; 988 size_t countI = 0;
982 double rc=0.0, rq=0.0; 989 double rc=0.0, rq=0.0;
983 do { 990 do {
984 double x = *data; 991 double x = *data;
985 double t,y; 992 double t,y;
986 if (!isnan(x)) { 993 if (!ISNAN(x)) {
987 countI++; 994 countI++;
988 //sum += x; [1] 995 //sum += x; [1]
989 y = x-rc; 996 y = x-rc;
990 t = sum+y; 997 t = sum+y;
991 rc= (t-sum)-y; 998 rc= (t-sum)-y;