The R Project SVN – R – Diff – /trunk/src/main/array.c



    int type;
    Rboolean NaRm, keepNA;
 
    checkArity(op, args);
    x = CAR(args); args = CDR(args);
    R_xlen_t n = asVecSize(CAR(args)); args = CDR(args);
    R_xlen_t p = asVecSize(CAR(args)); args = CDR(args);
    NaRm = asLogical(CAR(args));
    if (n == NA_INTEGER || n < 0)
	error(_("invalid '%s' argument"), "n");
    if (p == NA_INTEGER || p < 0)
	error(_("invalid '%s' argument"), "p");

	else
	    nthreads = 1; /* for now */
#pragma omp parallel for num_threads(nthreads) default(none) \
    firstprivate(x, ans, n, p, type, NaRm, keepNA, R_NaReal, R_NaInt, OP)
#endif
	for (R_xlen_t j = 0; j < p; j++) {
	    R_xlen_t  cnt = n, i;
	    LDOUBLE sum = 0.0;
	    switch (type) {
	    case REALSXP:
	    {
		double *rx = REAL(x) + (R_xlen_t)n*j;

	    rans = (LDOUBLE *) alloca(n * sizeof(LDOUBLE));
	    Memzero(rans, n);
	} else rans = Calloc(n, LDOUBLE);
	if (!keepNA && OP == 3) Cnt = Calloc(n, int);
 
	for (R_xlen_t j = 0; j < p; j++) {
	    LDOUBLE *ra = rans;
	    switch (type) {
	    case REALSXP:
	    {
		double *rx = REAL(x) + (R_xlen_t)n * j;
		if (keepNA)
		    for (R_xlen_t i = 0; i < n; i++) *ra++ += *rx++;
		else
		    for (R_xlen_t i = 0; i < n; i++, ra++, rx++)
			if (!ISNAN(*rx)) {
			    *ra += *rx;
			    if (OP == 3) Cnt[i]++;
			}
		break;
	    }
	    case INTSXP:
	    {
		int *ix = INTEGER(x) + (R_xlen_t)n * j;
		for (R_xlen_t i = 0; i < n; i++, ra++, ix++)
		    if (keepNA) {
			if (*ix != NA_INTEGER) *ra += *ix;
			else *ra = NA_REAL;
		    }
		    else if (*ix != NA_INTEGER) {

		break;
	    }
	    case LGLSXP:
	    {
		int *ix = LOGICAL(x) + (R_xlen_t)n * j;
		for (R_xlen_t i = 0; i < n; i++, ra++, ix++)
		    if (keepNA) {
			if (*ix != NA_LOGICAL) *ra += *ix;
			else *ra = NA_REAL;
		    }
		    else if (*ix != NA_LOGICAL) {

	    }
	    }
	}
	if (OP == 3) {
	    if (keepNA)
		for (R_xlen_t i = 0; i < n; i++) rans[i] /= p;
	    else
		for (R_xlen_t i = 0; i < n; i++) rans[i] /= Cnt[i];
	}
	for (R_xlen_t i = 0; i < n; i++) REAL(ans)[i] = (double) rans[i];
 
	if (!keepNA && OP == 3) Free(Cnt);
	if(n > 10000) Free(rans);
    }
 

Rev 64246	Rev 64374
Line 1254...	Line 1254...
1254	`int type;`	1254	`int type;`
1255	`Rboolean NaRm, keepNA;`	1255	`Rboolean NaRm, keepNA;`
1256		1256
1257	`checkArity(op, args);`	1257	`checkArity(op, args);`
1258	`x = CAR(args); args = CDR(args);`	1258	`x = CAR(args); args = CDR(args);`
1259	`int n = asInteger(CAR(args)); args = CDR(args);`	1259	`R_xlen_t n = asVecSize(CAR(args)); args = CDR(args);`
1260	`int p = asInteger(CAR(args)); args = CDR(args);`	1260	`R_xlen_t p = asVecSize(CAR(args)); args = CDR(args);`
1261	`NaRm = asLogical(CAR(args));`	1261	`NaRm = asLogical(CAR(args));`
1262	`if (n == NA_INTEGER \|\| n < 0)`	1262	`if (n == NA_INTEGER \|\| n < 0)`
1263	`error(_("invalid '%s' argument"), "n");`	1263	`error(_("invalid '%s' argument"), "n");`
1264	`if (p == NA_INTEGER \|\| p < 0)`	1264	`if (p == NA_INTEGER \|\| p < 0)`
1265	`error(_("invalid '%s' argument"), "p");`	1265	`error(_("invalid '%s' argument"), "p");`
Line 1285...	Line 1285...
1285	`else`	1285	`else`
1286	`nthreads = 1; /* for now */`	1286	`nthreads = 1; /* for now */`
1287	`#pragma omp parallel for num_threads(nthreads) default(none) \`	1287	`#pragma omp parallel for num_threads(nthreads) default(none) \`
1288	`firstprivate(x, ans, n, p, type, NaRm, keepNA, R_NaReal, R_NaInt, OP)`	1288	`firstprivate(x, ans, n, p, type, NaRm, keepNA, R_NaReal, R_NaInt, OP)`
1289	`#endif`	1289	`#endif`
1290	`for (int j = 0; j < p; j++) {`	1290	`for (R_xlen_t j = 0; j < p; j++) {`
1291	`int cnt = n, i;`	1291	`R_xlen_t cnt = n, i;`
1292	`LDOUBLE sum = 0.0;`	1292	`LDOUBLE sum = 0.0;`
1293	`switch (type) {`	1293	`switch (type) {`
1294	`case REALSXP:`	1294	`case REALSXP:`
1295	`{`	1295	`{`
1296	`double rx = REAL(x) + (R_xlen_t)nj;`	1296	`double rx = REAL(x) + (R_xlen_t)nj;`
Line 1336...	Line 1336...
1336	`rans = (LDOUBLE ) alloca(n sizeof(LDOUBLE));`	1336	`rans = (LDOUBLE ) alloca(n sizeof(LDOUBLE));`
1337	`Memzero(rans, n);`	1337	`Memzero(rans, n);`
1338	`} else rans = Calloc(n, LDOUBLE);`	1338	`} else rans = Calloc(n, LDOUBLE);`
1339	`if (!keepNA && OP == 3) Cnt = Calloc(n, int);`	1339	`if (!keepNA && OP == 3) Cnt = Calloc(n, int);`
1340		1340
1341	`for (int j = 0; j < p; j++) {`	1341	`for (R_xlen_t j = 0; j < p; j++) {`
1342	`LDOUBLE *ra = rans;`	1342	`LDOUBLE *ra = rans;`
1343	`switch (type) {`	1343	`switch (type) {`
1344	`case REALSXP:`	1344	`case REALSXP:`
1345	`{`	1345	`{`
1346	`double rx = REAL(x) + (R_xlen_t)n j;`	1346	`double rx = REAL(x) + (R_xlen_t)n j;`
1347	`if (keepNA)`	1347	`if (keepNA)`
1348	`for (int i = 0; i < n; i++) ra++ += rx++;`	1348	`for (R_xlen_t i = 0; i < n; i++) ra++ += rx++;`
1349	`else`	1349	`else`
1350	`for (int i = 0; i < n; i++, ra++, rx++)`	1350	`for (R_xlen_t i = 0; i < n; i++, ra++, rx++)`
1351	`if (!ISNAN(*rx)) {`	1351	`if (!ISNAN(*rx)) {`
1352	`ra += rx;`	1352	`ra += rx;`
1353	`if (OP == 3) Cnt[i]++;`	1353	`if (OP == 3) Cnt[i]++;`
1354	`}`	1354	`}`
1355	`break;`	1355	`break;`
1356	`}`	1356	`}`
1357	`case INTSXP:`	1357	`case INTSXP:`
1358	`{`	1358	`{`
1359	`int ix = INTEGER(x) + (R_xlen_t)n j;`	1359	`int ix = INTEGER(x) + (R_xlen_t)n j;`
1360	`for (int i = 0; i < n; i++, ra++, ix++)`	1360	`for (R_xlen_t i = 0; i < n; i++, ra++, ix++)`
1361	`if (keepNA) {`	1361	`if (keepNA) {`
1362	`if (ix != NA_INTEGER) ra += *ix;`	1362	`if (ix != NA_INTEGER) ra += *ix;`
1363	`else *ra = NA_REAL;`	1363	`else *ra = NA_REAL;`
1364	`}`	1364	`}`
1365	`else if (*ix != NA_INTEGER) {`	1365	`else if (*ix != NA_INTEGER) {`
Line 1369...	Line 1369...
1369	`break;`	1369	`break;`
1370	`}`	1370	`}`
1371	`case LGLSXP:`	1371	`case LGLSXP:`
1372	`{`	1372	`{`
1373	`int ix = LOGICAL(x) + (R_xlen_t)n j;`	1373	`int ix = LOGICAL(x) + (R_xlen_t)n j;`
1374	`for (int i = 0; i < n; i++, ra++, ix++)`	1374	`for (R_xlen_t i = 0; i < n; i++, ra++, ix++)`
1375	`if (keepNA) {`	1375	`if (keepNA) {`
1376	`if (ix != NA_LOGICAL) ra += *ix;`	1376	`if (ix != NA_LOGICAL) ra += *ix;`
1377	`else *ra = NA_REAL;`	1377	`else *ra = NA_REAL;`
1378	`}`	1378	`}`
1379	`else if (*ix != NA_LOGICAL) {`	1379	`else if (*ix != NA_LOGICAL) {`
Line 1384...	Line 1384...
1384	`}`	1384	`}`
1385	`}`	1385	`}`
1386	`}`	1386	`}`
1387	`if (OP == 3) {`	1387	`if (OP == 3) {`
1388	`if (keepNA)`	1388	`if (keepNA)`
1389	`for (int i = 0; i < n; i++) rans[i] /= p;`	1389	`for (R_xlen_t i = 0; i < n; i++) rans[i] /= p;`
1390	`else`	1390	`else`
1391	`for (int i = 0; i < n; i++) rans[i] /= Cnt[i];`	1391	`for (R_xlen_t i = 0; i < n; i++) rans[i] /= Cnt[i];`
1392	`}`	1392	`}`
1393	`for (int i = 0; i < n; i++) REAL(ans)[i] = (double) rans[i];`	1393	`for (R_xlen_t i = 0; i < n; i++) REAL(ans)[i] = (double) rans[i];`
1394		1394
1395	`if (!keepNA && OP == 3) Free(Cnt);`	1395	`if (!keepNA && OP == 3) Free(Cnt);`
1396	`if(n > 10000) Free(rans);`	1396	`if(n > 10000) Free(rans);`
1397	`}`	1397	`}`
1398		1398

The R Project SVN R

(root)/trunk/src/main/array.c @ 63243 – Rev 64246 → 64374