The R Project SVN R

Rev

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

Rev 36990 Rev 37393
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  R : A Computer Language for Statistical Data Analysis
 2 
 *  R : A Computer Language for Statistical Data Analysis
3 
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
 3 
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
4 
 *  Copyright (C) 1998-2001   The R Development Core Team
 4 
 *  Copyright (C) 1998-2001   The R Development Core Team
5 
 *  Copyright (C) 2002--2005  The R Foundation
 5 
 *  Copyright (C) 2002--2006  The R Foundation
6 
 *
 6 
 *
7 
 *  This program is free software; you can redistribute it and/or modify
 7 
 *  This program is free software; you can redistribute it and/or modify
8 
 *  it under the terms of the GNU General Public License as published by
 8 
 *  it under the terms of the GNU General Public License as published by
9 
 *  the Free Software Foundation; either version 2 of the License, or
 9 
 *  the Free Software Foundation; either version 2 of the License, or
10 
 *  (at your option) any later version.
 10 
 *  (at your option) any later version.
Line 386... Line 386...
386 
static void matprod(double *x, int nrx, int ncx,
 386 
static void matprod(double *x, int nrx, int ncx,
387 
		    double *y, int nry, int ncy, double *z)
 387 
		    double *y, int nry, int ncy, double *z)
388 
{
 388 
{
389 
    char *transa = "N", *transb = "N";
 389 
    char *transa = "N", *transb = "N";
390 
    int i,  j, k;
 390 
    int i,  j, k;
391 
    double one = 1.0, zero = 0.0, sum;
 391 
    double one = 1.0, zero = 0.0;
- 
 
 392 
    LDOUBLE sum;
392 
    Rboolean have_na = FALSE;
 393 
    Rboolean have_na = FALSE;
393 
 
 394 
 
394 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
 395 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
395 
	/* Don't trust the BLAS to handle NA/NaNs correctly: PR#4582
 396 
	/* Don't trust the BLAS to handle NA/NaNs correctly: PR#4582
396 
	 * The test is only O(n) here
 397 
	 * The test is only O(n) here
Line 430... Line 431...
430 
    } else { /* zero-extent operations should return zeroes */
 431 
    } else { /* zero-extent operations should return zeroes */
431 
	for(i = 0; i < nrx*ncy; i++) z[i].r = z[i].i = 0;
 432 
	for(i = 0; i < nrx*ncy; i++) z[i].r = z[i].i = 0;
432 
    }
 433 
    }
433 
#else
 434 
#else
434 
    int i, j, k;
 435 
    int i, j, k;
435 
    double xij_r, xij_i, yjk_r, yjk_i, sum_i, sum_r;
 436 
    double xij_r, xij_i, yjk_r, yjk_i;
- 
 
 437 
    LDOUBLE sum_i, sum_r;
436 
 
 438 
 
437 
    for (i = 0; i < nrx; i++)
 439 
    for (i = 0; i < nrx; i++)
438 
	for (k = 0; k < ncy; k++) {
 440 
	for (k = 0; k < ncy; k++) {
439 
	    z[i + k * nrx].r = NA_REAL;
 441 
	    z[i + k * nrx].r = NA_REAL;
440 
	    z[i + k * nrx].i = NA_REAL;
 442 
	    z[i + k * nrx].i = NA_REAL;
Line 1115... Line 1117...
1115 
{
 1117 
{
1116 
    SEXP x, ans = R_NilValue;
 1118 
    SEXP x, ans = R_NilValue;
1117 
    int OP, n, p, cnt = 0, i, j, type;
 1119 
    int OP, n, p, cnt = 0, i, j, type;
1118 
    Rboolean NaRm, keepNA;
 1120 
    Rboolean NaRm, keepNA;
1119 
    int *ix;
 1121 
    int *ix;
- 
 
 1122 
    double *rx;
1120 
    double *rx, sum = 0.0;
 1123 
    LDOUBLE sum = 0.0;
1121 
 
 1124 
 
1122 
    checkArity(op, args);
 1125 
    checkArity(op, args);
1123 
    x = CAR(args); args = CDR(args);
 1126 
    x = CAR(args); args = CDR(args);
1124 
    n = asInteger(CAR(args)); args = CDR(args);
 1127 
    n = asInteger(CAR(args)); args = CDR(args);
1125 
    p = asInteger(CAR(args)); args = CDR(args);
 1128 
    p = asInteger(CAR(args)); args = CDR(args);
Line 1180... Line 1183...
1180 
 
 1183 
 
1181 
    if (OP == 2 || OP == 3) { /* rows */
 1184 
    if (OP == 2 || OP == 3) { /* rows */
1182 
	cnt = p;
 1185 
	cnt = p;
1183 
	PROTECT(ans = allocVector(REALSXP, n));
 1186 
	PROTECT(ans = allocVector(REALSXP, n));
1184 
 
 1187 
 
1185 
	/* reverse summation order to improve cache hits */
 1188 
	/* interchange summation order to improve cache hits */
1186 
	if (type == REALSXP) {
 1189 
	if (type == REALSXP) {
1187 
	    double *rans = REAL(ans), *ra = rans;
 - 
 
1188 
	    int *Cnt = NULL, *c;
 1190 
	    int *Cnt = NULL, *c;
- 
 
 1191 
	    LDOUBLE *rans, *ra;
- 
 
 1192 
	    if(n <= 10000) {
- 
 
 1193 
		rans = (LDOUBLE *) alloca(n * sizeof(LDOUBLE));
- 
 
 1194 
		R_CheckStack();
- 
 
 1195 
		memset(rans, 0, n*sizeof(LDOUBLE));
- 
 
 1196 
	    } else rans = Calloc(n, LDOUBLE);
1189 
	    rx = REAL(x);
 1197 
	    rx = REAL(x);
1190 
	    if (!keepNA && OP == 3) Cnt = Calloc(n, int);
 1198 
	    if (!keepNA && OP == 3) Cnt = Calloc(n, int);
1191 
	    memset(rans, 0, n*sizeof(double));
 - 
 
1192 
	    for (j = 0; j < p; j++) {
 1199 
	    for (j = 0; j < p; j++) {
1193 
		ra = rans;
 1200 
		ra = rans;
1194 
		if (keepNA)
 1201 
		if (keepNA)
1195 
		    for (i = 0; i < n; i++) *ra++ += *rx++;
 1202 
		    for (i = 0; i < n; i++) *ra++ += *rx++;
1196 
		else
 1203 
		else
Line 1207... Line 1214...
1207 
		    for (ra = rans, c = Cnt, i = 0; i < n; i++, c++)
 1214 
		    for (ra = rans, c = Cnt, i = 0; i < n; i++, c++)
1208 
			if (*c > 0) *ra++ /= *c; else *ra++ = NA_REAL;
 1215 
			if (*c > 0) *ra++ /= *c; else *ra++ = NA_REAL;
1209 
		    Free(Cnt);
 1216 
		    Free(Cnt);
1210 
		}
 1217 
		}
1211 
	    }
 1218 
	    }
- 
 
 1219 
	    for (i = 0; i < n; i++) REAL(ans)[i] = rans[i];
- 
 
 1220 
	    if(n > 10000) Free(rans);
1212 
	    UNPROTECT(1);
 1221 
	    UNPROTECT(1);
1213 
	    return ans;
 1222 
	    return ans;
1214 
	}
 1223 
	}
1215 
 
 1224 
 
1216 
	for (i = 0; i < n; i++) {
 1225 
	for (i = 0; i < n; i++) {
1217 
	    switch (type) {
 1226 
	    switch (type) {
1218 
	    case REALSXP: /* this cannot be reached */
 - 
 
1219 
		rx = REAL(x) + i;
 - 
 
1220 
		if (keepNA)
 - 
 
1221 
		    for (sum = 0., j = 0; j < p; j++, rx += n) sum += *rx;
 - 
 
1222 
		else {
 - 
 
1223 
		    for (cnt = 0, sum = 0., j = 0; j < p; j++, rx += n)
 - 
 
1224 
			if (!ISNAN(*rx)) {cnt++; sum += *rx;}
 - 
 
1225 
			else if (keepNA) {sum = NA_REAL; break;}
 - 
 
1226 
		}
 - 
 
1227 
		break;
 - 
 
1228 
	    case INTSXP:
 1227 
	    case INTSXP:
1229 
		ix = INTEGER(x) + i;
 1228 
		ix = INTEGER(x) + i;
1230 
		for (cnt = 0, sum = 0., j = 0; j < p; j++, ix += n)
 1229 
		for (cnt = 0, sum = 0., j = 0; j < p; j++, ix += n)
1231 
		    if (*ix != NA_INTEGER) {cnt++; sum += *ix;}
 1230 
		    if (*ix != NA_INTEGER) {cnt++; sum += *ix;}
1232 
		    else if (keepNA) {sum = NA_REAL; break;}
 1231 
		    else if (keepNA) {sum = NA_REAL; break;}