The R Project SVN R

Rev

Rev 17515 | Rev 18753 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 17515 Rev 17521
Line 932... Line 932...
932 
/* colSums(x, n, p, na.rm) and friends */
 932 
/* colSums(x, n, p, na.rm) and friends */
933 
SEXP do_colsum(SEXP call, SEXP op, SEXP args, SEXP rho)
 933 
SEXP do_colsum(SEXP call, SEXP op, SEXP args, SEXP rho)
934 
{
 934 
{
935 
    SEXP x, ans = R_NilValue;
 935 
    SEXP x, ans = R_NilValue;
936 
    int OP, n, p, cnt = 0, i, j, type;
 936 
    int OP, n, p, cnt = 0, i, j, type;
937 
    Rboolean NaRm;
 937 
    Rboolean NaRm, keepNA;
938 
    int *ix;
 938 
    int *ix;
939 
    double *rx, sum = 0.0;
 939 
    double *rx, sum = 0.0;
940 
    /* Rcomplex cx, csum;*/
 - 
 
941 
 
 940 
 
942 
    checkArity(op, args);
 941 
    checkArity(op, args);
943 
    x = CAR(args); args = CDR(args);
 942 
    x = CAR(args); args = CDR(args);
944 
    n = asInteger(CAR(args)); args = CDR(args);
 943 
    n = asInteger(CAR(args)); args = CDR(args);
945 
    p = asInteger(CAR(args)); args = CDR(args);
 944 
    p = asInteger(CAR(args)); args = CDR(args);
946 
    NaRm = asLogical(CAR(args));
 945 
    NaRm = asLogical(CAR(args));
947 
    if(n == NA_INTEGER || n <= 0)
 946 
    if (n == NA_INTEGER || n <= 0)
948 
	errorcall(call, "invalid value of n");
 947 
	errorcall(call, "invalid value of n");
949 
    if(p == NA_INTEGER || p <= 0)
 948 
    if (p == NA_INTEGER || p <= 0)
950 
	errorcall(call, "invalid value of p");
 949 
	errorcall(call, "invalid value of p");
951 
    if(NaRm == NA_LOGICAL) errorcall(call, "invalid value of na.rm");
 950 
    if (NaRm == NA_LOGICAL) errorcall(call, "invalid value of na.rm");
- 
 
 951 
    keepNA = !NaRm;
952 
 
 952 
 
953 
    OP = PRIMVAL(op);
 953 
    OP = PRIMVAL(op);
954 
    switch (type = TYPEOF(x)) {
 954 
    switch (type = TYPEOF(x)) {
955 
    case LGLSXP: break;
 955 
    case LGLSXP: break;
956 
    case INTSXP: break;
 956 
    case INTSXP: break;
957 
    case REALSXP: break;
 957 
    case REALSXP: break;
958 
/*    case CPLXSXP: break; */
 - 
 
959 
    default:
 958 
    default:
960 
	errorcall(call, "`x' must be numeric");
 959 
	errorcall(call, "`x' must be numeric");
961 
    }
 960 
    }
962 
 
 961 
 
963 
    if(OP == 0 || OP == 1) { /* columns */
 962 
    if (OP == 0 || OP == 1) { /* columns */
- 
 
 963 
	cnt = n;
964 
	PROTECT(ans = allocVector(REALSXP, p));
 964 
	PROTECT(ans = allocVector(REALSXP, p));
965 
	for (j = 0; j < p; j++) {
 965 
	for (j = 0; j < p; j++) {
966 
	    switch (type) {
 966 
	    switch (type) {
967 
	    case REALSXP:
 967 
	    case REALSXP:
968 
		rx = REAL(x);
 968 
		rx = REAL(x) + n*j;
- 
 
 969 
#ifdef IEEE_754
- 
 
 970 
		if (keepNA)
- 
 
 971 
		    for (sum = 0., i = 0; i < n; i++) sum += *rx++;
- 
 
 972 
		else {
- 
 
 973 
		    for (cnt = 0, sum = 0., i = 0; i < n; i++, rx++)
- 
 
 974 
			if (!ISNAN(*rx)) {cnt++; sum += *rx;}
- 
 
 975 
			else if (keepNA) {sum = NA_REAL; break;}
- 
 
 976 
		}
- 
 
 977 
#else
969 
		for (cnt = 0, sum = 0., i = 0; i < n; i++)
 978 
		for (cnt = 0, sum = 0., i = 0; i < n; i++, rx++)
970 
		    if(!ISNAN(rx[i+n*j])) {cnt++; sum += rx[i+n*j];}
 979 
		    if (!ISNAN(*rx)) {cnt++; sum += *rx;}
971 
		    else if(!NaRm) {sum = NA_REAL; break;}
 980 
		    else if (keepNA) {sum = NA_REAL; break;}
- 
 
 981 
#endif
972 
		break;
 982 
		break;
973 
	    case INTSXP:
 983 
	    case INTSXP:
974 
		ix = INTEGER(x);
 984 
		ix = INTEGER(x) + n*j;
975 
		for (cnt = 0, sum = 0., i = 0; i < n; i++)
 985 
		for (cnt = 0, sum = 0., i = 0; i < n; i++, ix++)
976 
		    if(ix[i+n*j] != NA_INTEGER) {cnt++; sum += ix[i+n*j];}
 986 
		    if (*ix != NA_INTEGER) {cnt++; sum += *ix;}
977 
		    else if(!NaRm) {sum = NA_REAL; break;}
 987 
		    else if (keepNA) {sum = NA_REAL; break;}
978 
		break;
 988 
		break;
979 
	    case LGLSXP:
 989 
	    case LGLSXP:
980 
		ix = LOGICAL(x);
 990 
		ix = LOGICAL(x) + n*j;
981 
		for (cnt = 0, sum = 0., i = 0; i < n; i++)
 991 
		for (cnt = 0, sum = 0., i = 0; i < n; i++, ix++)
982 
		    if(ix[i+n*j] != NA_LOGICAL) {cnt++; sum += ix[i+n*j];}
 992 
		    if (*ix != NA_LOGICAL) {cnt++; sum += *ix;}
983 
		    else if(!NaRm) {sum = NA_REAL; break;}
 993 
		    else if (keepNA) {sum = NA_REAL; break;}
984 
		break;
 994 
		break;
985 
	    }
 995 
	    }
986 
	    if(OP == 1) {
 996 
	    if (OP == 1) {
987 
		if(cnt > 0) sum /= cnt; else sum = NA_REAL;
 997 
		if (cnt > 0) sum /= cnt; else sum = NA_REAL;
988 
	    }
 998 
	    }
989 
	    REAL(ans)[j] = sum;
 999 
	    REAL(ans)[j] = sum;
990 
	}
 1000 
	}
991 
    }
 1001 
    }
992 
 
 1002 
 
993 
    if(OP == 2 || OP == 3) { /* rows */
 1003 
    if (OP == 2 || OP == 3) { /* rows */
- 
 
 1004 
	cnt = p;
994 
	PROTECT(ans = allocVector(REALSXP, n));
 1005 
	PROTECT(ans = allocVector(REALSXP, n));
- 
 
 1006 
 
- 
 
 1007 
#ifdef IEEE_754
- 
 
 1008 
	/* reverse summation order to improve cache hits */
- 
 
 1009 
	if (type == REALSXP) {
- 
 
 1010 
	    double *rans = REAL(ans), *ra = rans, *cnt = NULL, *c;
- 
 
 1011 
	    rx = REAL(x);
- 
 
 1012 
	    if (!keepNA && OP == 3) cnt = Calloc(n, double);
- 
 
 1013 
	    for (ra = rans, i = 0; i < n; i++) *ra++ = 0.0;
- 
 
 1014 
	    for (j = 0; j < p; j++) {
- 
 
 1015 
		ra = rans;
- 
 
 1016 
		if (keepNA)
- 
 
 1017 
		    for (i = 0; i < n; i++) *ra++ += *rx++;
- 
 
 1018 
		else
- 
 
 1019 
		    for (c = cnt, i = 0; i < n; i++, ra++, rx++, c++)
- 
 
 1020 
			if (!ISNAN(*rx)) {
- 
 
 1021 
			    *ra += *rx;
- 
 
 1022 
			    if (OP == 3) (*c)++;
- 
 
 1023 
			}
- 
 
 1024 
	    }
- 
 
 1025 
	    if (OP == 3) {
- 
 
 1026 
		if (keepNA)
- 
 
 1027 
		    for (ra = rans, i = 0; i < n; i++)
- 
 
 1028 
			*ra++ /= p;
- 
 
 1029 
		else {
- 
 
 1030 
		    for (ra = rans, c = cnt, i = 0; i < n; i++, c++)
- 
 
 1031 
			if (*c > 0) *ra++ /= *c; else *ra++ = NA_REAL;
- 
 
 1032 
		    Free(cnt);
- 
 
 1033 
		}
- 
 
 1034 
	    }
- 
 
 1035 
	    UNPROTECT(1);
- 
 
 1036 
	    return ans;
- 
 
 1037 
	}
- 
 
 1038 
#endif
- 
 
 1039 
 
995 
	for (i = 0; i < n; i++) {
 1040 
	for (i = 0; i < n; i++) {
996 
	    switch (type) {
 1041 
	    switch (type) {
997 
	    case REALSXP:
 1042 
	    case REALSXP:
998 
		rx = REAL(x);
 1043 
		rx = REAL(x) + i;
- 
 
 1044 
#ifdef IEEE_754
- 
 
 1045 
		if (keepNA)
- 
 
 1046 
		    for (sum = 0., j = 0; j < p; j++, rx += n) sum += *rx;
- 
 
 1047 
		else {
- 
 
 1048 
		    for (cnt = 0, sum = 0., j = 0; j < p; j++, rx += n)
- 
 
 1049 
			if (!ISNAN(*rx)) {cnt++; sum += *rx;}
- 
 
 1050 
			else if (keepNA) {sum = NA_REAL; break;}
- 
 
 1051 
		}
- 
 
 1052 
#else
999 
		for (cnt = 0, sum = 0., j = 0; j < p; j++)
 1053 
		for (cnt = 0, sum = 0., j = 0; j < p; j++, rx += n)
1000 
		    if(!ISNAN(rx[i+n*j])) {cnt++; sum += rx[i+n*j];}
 1054 
		    if (!ISNAN(*rx)) {cnt++; sum += *rx;}
1001 
		    else if(!NaRm) {sum = NA_REAL; break;}
 1055 
		    else if (keepNA) {sum = NA_REAL; break;}
- 
 
 1056 
#endif
1002 
		break;
 1057 
		break;
1003 
	    case INTSXP:
 1058 
	    case INTSXP:
1004 
		ix = INTEGER(x);
 1059 
		ix = INTEGER(x) + i;
1005 
		for (cnt = 0, sum = 0., j = 0; j < p; j++)
 1060 
		for (cnt = 0, sum = 0., j = 0; j < p; j++, ix += n)
1006 
		    if(ix[i+n*j] != NA_INTEGER) {cnt++; sum += ix[i+n*j];}
 1061 
		    if (*ix != NA_INTEGER) {cnt++; sum += *ix;}
1007 
		    else if(!NaRm) {sum = NA_REAL; break;}
 1062 
		    else if (keepNA) {sum = NA_REAL; break;}
1008 
		break;
 1063 
		break;
1009 
	    case LGLSXP:
 1064 
	    case LGLSXP:
1010 
		ix = LOGICAL(x);
 1065 
		ix = LOGICAL(x) + i;
1011 
		for (cnt = 0, sum = 0., j = 0; j < p; j++)
 1066 
		for (cnt = 0, sum = 0., j = 0; j < p; j++, ix += n)
1012 
		    if(ix[i+n*j] != NA_LOGICAL) {cnt++; sum += ix[i+n*j];}
 1067 
		    if (*ix != NA_LOGICAL) {cnt++; sum += *ix;}
1013 
		    else if(!NaRm) {sum = NA_REAL; break;}
 1068 
		    else if (keepNA) {sum = NA_REAL; break;}
1014 
		break;
 1069 
		break;
1015 
	    }
 1070 
	    }
1016 
	    if(OP == 3) {
 1071 
	    if (OP == 3) {
1017 
		if(cnt > 0) sum /= cnt; else sum = NA_REAL;
 1072 
		if (cnt > 0) sum /= cnt; else sum = NA_REAL;
1018 
	    }
 1073 
	    }
1019 
	    REAL(ans)[i] = sum;
 1074 
	    REAL(ans)[i] = sum;
1020 
	}
 1075 
	}
1021 
    }
 1076 
    }
1022 
 
 1077