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

Rev	Author	Line No.	Line
2	r	1	`/*`
		2	`* R : A Computer Langage for Statistical Data Analysis`
		3	`* Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka`
		4	`*`
		5	`* This program is free software; you can redistribute it and/or modify`
		6	`* it under the terms of the GNU General Public License as published by`
		7	`* the Free Software Foundation; either version 2 of the License, or`
		8	`* (at your option) any later version.`
		9	`*`
		10	`* This program is distributed in the hope that it will be useful,`
		11	`* but WITHOUT ANY WARRANTY; without even the implied warranty of`
		12	`* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the`
		13	`* GNU General Public License for more details.`
		14	`*`
		15	`* You should have received a copy of the GNU General Public License`
		16	`* along with this program; if not, write to the Free Software`
		17	`* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.`
		18	`*/`
		19
		20	`#include "Defn.h"`
		21	`#include "Mathlib.h"`
		22
		23	`static void CheckDims(SEXP dims)`
		24	`{`
		25	`int i;`
		26
		27	`for (i = 0; i < LENGTH(dims); i++) {`
		28	`if (INTEGER(dims)[i] <= 0)`
		29	`error("invalid array extent\n");`
		30	`}`
		31	`}`
		32
		33	`SEXP do_matrix(SEXP call, SEXP op, SEXP args, SEXP rho)`
		34	`{`
		35	`SEXP vals, snr, snc;`
		36	`int nr, nc, byrow, lendat;`
		37
		38	`checkArity(op, args);`
		39	`vals = CAR(args);`
		40	`snr = CADR(args);`
		41	`snc = CADDR(args);`
		42	`byrow = asInteger(CADR(CDDR(args)));`
		43
		44	`if (isVector(vals) \|\| isList(vals)) {`
		45	`if(length(vals) <= 0)`
		46	`errorcall(call, "argument has length zero\n");`
		47	`} else errorcall(call, "invalid matrix element type\n");`
		48
		49	`if (!isNumeric(snr) \|\| !isNumeric(snc))`
		50	`error("non-numeric matrix extent\n");`
		51
		52	`lendat = length(vals);`
		53	`nr = asInteger(snr);`
		54	`nc = asInteger(snc);`
		55
		56
		57	`if( lendat != 1 && (nr*nc) % lendat != 0 ) {`
		58	`if( ((lendat>nr) && (lendat/nr)*nr != lendat ) \|\|`
		59	`((lendat< nr) && (nr/lendat) * lendat != nr ))`
		60	`warning("Replacement length not a multiple of the elements to replace in matrix(...) \n");`
		61	`else if( ((lendat>nc) && (lendat/nc)*nc != lendat ) \|\|`
		62	`((lendat< nc) && (nc/lendat) * lendat != nc ))`
		63	`warning("Replacement length not a multiple of the elements to replace in matrix(...) \n");`
		64	`}`
		65
		66	`PROTECT(snr = allocMatrix(TYPEOF(vals), nr, nc));`
		67	`LEVELS(snr) = LEVELS(vals);`
		68	`if(isVector(vals))`
		69	`copyMatrix(snr, vals, byrow);`
		70	`else`
		71	`copyListMatrix(snr, vals, byrow);`
		72	`UNPROTECT(1);`
		73	`return snr;`
		74	`}`
		75
		76
		77	`SEXP allocMatrix(SEXPTYPE mode, int nrow, int ncol)`
		78	`{`
		79	`SEXP s, t;`
		80	`int n;`
		81
		82	`if (nrow <= 0 \|\| ncol <= 0)`
		83	`error("nonpositive extents to matrix\n");`
		84	`n = nrow * ncol;`
		85	`PROTECT(s = allocVector(mode, n));`
		86	`PROTECT(t = allocVector(INTSXP, 2));`
		87	`INTEGER(t)[0] = nrow;`
		88	`INTEGER(t)[1] = ncol;`
		89	`setAttrib(s, R_DimSymbol, t);`
		90	`UNPROTECT(2);`
		91	`return s;`
		92	`}`
		93
		94	`SEXP do_array(SEXP call, SEXP op, SEXP args, SEXP rho)`
		95	`{`
		96	`SEXP vals, dims, ans;`
		97
		98	`if (length(args) != 2)`
		99	`error("incorrect arg count to \"array\"\n");`
		100	`vals = CAR(args);`
		101	`dims = CADR(args);`
		102
		103	`if (isVector(vals) && isNumeric(dims) && LENGTH(dims) >= 2) {`
		104	`PROTECT(dims = coerceVector(dims, INTSXP));`
		105	`CheckDims(dims);`
		106	`PROTECT(ans = allocArray(TYPEOF(vals), dims));`
		107	`LEVELS(ans) = LEVELS(vals);`
		108	`copyVector(ans, vals);`
		109	`UNPROTECT(2);`
		110	`return ans;`
		111	`}`
		112	`else`
		113	`error("bad arguments to array\n");`
		114	`/NOTREACHED/`
		115	`}`
		116
		117	`SEXP allocArray(SEXPTYPE mode, SEXP dims)`
		118	`{`
		119	`SEXP array;`
		120	`int i, n;`
		121
		122	`n = 1;`
		123	`for (i = 0; i < LENGTH(dims); i++)`
		124	`n = n * INTEGER(dims)[i];`
		125
		126	`PROTECT(dims = duplicate(dims));`
		127	`PROTECT(array = allocVector(mode, n));`
		128	`setAttrib(array, R_DimSymbol, dims);`
		129	`UNPROTECT(2);`
		130	`return array;`
		131	`}`
		132
		133	`/* DropDims strips away redundant dimensioning */`
		134	`/* information If there is an appropriate dimnames */`
		135	`/* attribute the correct element is extracted and */`
		136	`/* attached to the vector as a names attribute. Note */`
		137	`/* that this function mutates x. Duplication should */`
		138	`/* occur before this is called. */`
		139
		140	`SEXP DropDims(SEXP x)`
		141	`{`
		142	`SEXP p, q, r, pdims, pdimnames;`
		143	`int i, n;`
		144
		145	`/* Locate the dims and dimnames attributes */`
		146
		147	`pdims = R_NilValue;`
		148	`pdimnames = R_NilValue;`
		149	`for(p=ATTRIB(x) ; p!=R_NilValue ; p=CDR(p)) {`
		150	`if(TAG(p) == R_DimSymbol)`
		151	`pdims = p;`
		152	`if(TAG(p) == R_DimNamesSymbol)`
		153	`pdimnames = p;`
		154	`}`
		155
		156	`/* Check that dropping will actually do something */`
		157
		158	`if(pdims == R_NilValue) return x;`
		159	`p = CAR(pdims);`
		160	`n = 0;`
		161	`for(i=0 ; i<LENGTH(p) ; i++)`
		162	`if(INTEGER(p)[i] != 1) n++;`
		163	`if(n == LENGTH(p)) return x;`
		164
		165	`PROTECT(x);`
		166
		167	`if(n <= 1) { /* vector */`
		168	`SEXP newnames = R_NilValue;`
		169	`if(pdimnames != R_NilValue) {`
		170	`q = CAR(pdimnames);`
		171	`for(i=0 ; i<LENGTH(p) ; i++) {`
		172	`if(INTEGER(p)[i] != 1)`
		173	`newnames = CAR(q);`
		174	`q = CDR(q);`
		175	`}`
		176	`}`
		177	`PROTECT(newnames);`
		178	`setAttrib(x, R_DimNamesSymbol, R_NilValue);`
		179	`setAttrib(x, R_DimSymbol, R_NilValue);`
		180	`setAttrib(x, R_NamesSymbol, newnames);`
		181	`UNPROTECT(1);`
		182	`}`
		183	`else { /* array */`
		184	`SEXP newdims, newdimnames;`
		185	`PROTECT(newdims = allocVector(INTSXP, n));`
		186	`p = CAR(pdims);`
		187	`n = 0;`
		188	`for(i=0 ; i<LENGTH(p) ; i++)`
		189	`if(INTEGER(p)[i] != 1)`
		190	`INTEGER(newdims)[n++] = INTEGER(p)[i];`
		191	`if(pdimnames) {`
		192	`PROTECT(newdimnames = allocList(n));`
		193	`q = CAR(pdimnames);`
		194	`r = newdimnames;`
		195	`for(i=0 ; i<LENGTH(p) ; i++) {`
		196	`if(INTEGER(p)[i] != 1) {`
		197	`CAR(r) = CAR(q);`
		198	`r = CDR(r);`
		199	`}`
		200	`q = CDR(q);`
		201	`}`
		202	`UNPROTECT(1);`
		203	`}`
		204	`CAR(pdims) = newdims;`
		205	`if(pdimnames != R_NilValue)`
		206	`CAR(pdimnames) = newdimnames;`
		207	`UNPROTECT(1);`
		208	`}`
		209	`UNPROTECT(1);`
		210	`return x;`
		211	`}`
		212
		213	`SEXP do_drop(SEXP call, SEXP op, SEXP args, SEXP rho)`
		214	`{`
		215	`SEXP x, xdims;`
		216	`int i, n, shorten;`
		217
		218	`checkArity(op, args);`
		219	`x = CAR(args);`
		220	`if((xdims = getAttrib(x, R_DimSymbol)) != R_NilValue) {`
		221	`n = LENGTH(xdims);`
		222	`shorten = 0;`
		223	`for(i=0 ; i<n ; i++)`
		224	`if(INTEGER(xdims)[i] <= 1) shorten = 1;`
		225	`if(shorten) {`
		226	`if(NAMED(x)) x = duplicate(x);`
		227	`x = DropDims(x);`
		228	`}`
		229	`}`
		230	`return x;`
		231	`}`
		232
		233	`/* Length of Primitive Objects */`
		234
		235	`SEXP do_length(SEXP call, SEXP op, SEXP args, SEXP rho)`
		236	`{`
		237	`SEXP ans;`
		238
		239	`if (length(args) != 1)`
		240	`error("incorrect number of args to length\n");`
		241
		242	`ans = allocVector(INTSXP, 1);`
		243
		244	`#ifdef OLD`
		245	`switch(TYPEOF(CAR(args))) {`
		246	`case NILSXP:`
		247	`INTEGER(ans)[0] = 0;`
		248	`break;`
		249	`case LGLSXP:`
		250	`case FACTSXP:`
		251	`case ORDSXP:`
		252	`case INTSXP:`
		253	`case REALSXP:`
		254	`case CPLXSXP:`
		255	`case STRSXP:`
		256	`case EXPRSXP:`
		257	`INTEGER(ans)[0] = LENGTH(CAR(args));`
		258	`break;`
		259	`case LISTSXP:`
		260	`case LANGSXP:`
		261	`INTEGER(ans)[0] = length(CAR(args));`
		262	`break;`
		263	`case ENVSXP:`
		264	`INTEGER(ans)[0] = length(FRAME(CAR(args)));`
		265	`break;`
		266	`default:`
		267	`INTEGER(ans)[0] = 1;`
		268	`break;`
		269	`}`
		270	`#else`
		271	`INTEGER(ans)[0] = length(CAR(args));`
		272	`#endif`
		273	`return ans;`
		274	`}`
		275
		276	`SEXP do_nlevels(SEXP call, SEXP op, SEXP args, SEXP rho)`
		277	`{`
		278	`SEXP ans;`
		279
		280	`checkArity(op, args);`
		281	`ans = allocVector(INTSXP, 1);`
		282	`if (isFactor(CAR(args)))`
		283	`INTEGER(ans)[0] = LEVELS(CAR(args));`
		284	`else`
		285	`INTEGER(ans)[0] = NA_INTEGER;`
		286	`return ans;`
		287	`}`
		288
		289
		290	`SEXP do_rowscols(SEXP call, SEXP op, SEXP args, SEXP rho)`
		291	`{`
		292	`SEXP ans;`
		293	`int i, j, nr, nc;`
		294
		295	`if (length(args) != 1)`
		296	`error("incorrect number of args to row/col\n");`
		297	`if (!isMatrix(CAR(args)))`
		298	`error("a matrix is required as arg to row/col\n");`
		299
		300	`nr = nrows(CAR(args));`
		301	`nc = ncols(CAR(args));`
		302
		303	`ans = allocMatrix(INTSXP, nr, nc);`
		304
		305	`switch (PRIMVAL(op)) {`
		306	`case 1:`
		307	`for (i = 0; i < nr; i++)`
		308	`for (j = 0; j < nc; j++)`
		309	`INTEGER(ans)[i + j * nr] = i + 1;`
		310	`break;`
		311	`case 2:`
		312	`for (i = 0; i < nr; i++)`
		313	`for (j = 0; j < nc; j++)`
		314	`INTEGER(ans)[i + j * nr] = j + 1;`
		315	`break;`
		316	`}`
		317	`return ans;`
		318	`}`
		319
		320	`static void matprod(double x, int nrx, int ncx, double y, int nry, int ncy, double *z)`
		321	`{`
		322	`int i, j, k;`
		323	`double xij, yjk, sum;`
		324
		325	`for (i = 0; i < nrx; i++)`
		326	`for (k = 0; k < ncy; k++) {`
		327	`z[i + k * nrx] = NA_REAL;`
		328	`sum = 0.0;`
		329	`for (j = 0; j < ncx; j++) {`
		330	`xij = x[i + j * nrx];`
		331	`yjk = y[j + k * nry];`
		332	`if (!FINITE(xij) \|\| !FINITE(yjk))`
		333	`goto next_ik;`
		334	`sum += xij * yjk;`
		335	`}`
		336	`z[i + k * nrx] = sum;`
		337	`next_ik:`
		338	`;`
		339	`}`
		340	`}`
		341
		342	`static void cmatprod(complex *x, int nrx, int ncx,`
		343	`complex y, int nry, int ncy, complex z)`
		344	`{`
		345	`int i, j, k;`
		346	`double xij_r, xij_i, yjk_r, yjk_i, sum_i, sum_r;`
		347
		348	`for (i=0; i<nrx; i++)`
		349	`for (k=0; k<ncy; k++) {`
		350	`z[i+k*nrx].r = NA_REAL;`
		351	`z[i+k*nrx].i = NA_REAL;`
		352	`sum_r = 0.0;`
		353	`sum_i = 0.0;`
		354	`for (j=0; j<ncx; j++) {`
		355	`xij_r = x[i+j*nrx].r;`
		356	`xij_i = x[i+j*nrx].i;`
		357	`yjk_r = y[j+k*nry].r;`
		358	`yjk_i = y[j+k*nry].i;`
		359	`if (!FINITE(xij_r) \|\| !FINITE(xij_i)`
		360	`\|\| !FINITE(yjk_r) \|\| !FINITE(yjk_i))`
		361	`goto next_ik;`
		362	`sum_r += (xij_r * yjk_r - xij_i * yjk_i);`
		363	`sum_i += (xij_r * yjk_i + xij_i * yjk_r);`
		364	`}`
		365	`z[i+k*nrx].r = sum_r;`
		366	`z[i+k*nrx].i = sum_i;`
		367	`next_ik:`
		368	`;`
		369	`}`
		370	`}`
		371
		372	`static void crossprod(double x, int nrx, int ncx, double y, int nry, int ncy, double *z)`
		373	`{`
		374	`int i, j, k;`
		375	`double xji, yjk, sum;`
		376
		377	`for (i = 0; i < ncx; i++)`
		378	`for (k = 0; k < ncy; k++) {`
		379	`z[i + k * ncx] = NA_REAL;`
		380	`sum = 0.0;`
		381	`for (j = 0; j < nrx; j++) {`
		382	`xji = x[j + i * nrx];`
		383	`yjk = y[j + k * nry];`
		384	`if (!FINITE(xji) \|\| !FINITE(yjk))`
		385	`goto next_ik;`
		386	`sum += xji * yjk;`
		387	`}`
		388	`z[i + k * ncx] = sum;`
		389	`next_ik:`
		390	`;`
		391	`}`
		392	`}`
		393
		394	`#ifdef COMPLEX_DATA`
		395	`static void ccrossprod(complex x, int nrx, int ncx, complex y, int nry, int ncy, complex *z)`
		396	`{`
		397	`int i, j, k;`
		398	`double xji_r, xji_i, yjk_r, yjk_i, sum_r, sum_i;`
		399
		400	`for (i = 0; i < ncx; i++)`
		401	`for (k = 0; k < ncy; k++) {`
		402	`z[i + k * ncx].r = NA_REAL;`
		403	`z[i + k * ncx].i = NA_REAL;`
		404	`sum_r = 0.0;`
		405	`sum_i = 0.0;`
		406	`for (j = 0; j < nrx; j++) {`
		407	`xji_r = x[j + i * nrx].r;`
		408	`xji_i = x[j + i * nrx].i;`
		409	`yjk_r = y[j + k * nry].r;`
		410	`yjk_i = y[j + k * nry].i;`
		411	`if (!FINITE(xji_r) \|\| !FINITE(xji_i)`
		412	`\|\| !FINITE(yjk_r) \|\| !FINITE(yjk_i))`
		413	`goto next_ik;`
		414	`sum_r += (xji_r * yjk_r - xji_i * yjk_i);`
		415	`sum_i += (xji_r * yjk_i + xji_i * yjk_r);`
		416	`}`
		417	`z[i + k * ncx].r = sum_r;`
		418	`z[i + k * ncx].i = sum_i;`
		419	`next_ik:`
		420	`;`
		421	`}`
		422	`}`
		423	`#endif`
		424
		425	`SEXP do_matprod(SEXP call, SEXP op, SEXP args, SEXP rho)`
		426	`{`
		427	`int ldx, ldy, nrx, ncx, nry, ncy, mode;`
		428	`SEXP x, y, xdims, ydims, ans;`
		429
		430	`#ifdef COMPLEX_DATA`
		431	`if (!(isNumeric(CAR(args)) \|\| isComplex(CAR(args))) \|\|`
		432	`!(isNumeric(CADR(args)) \|\| isComplex(CADR(args))))`
		433	`#else`
		434	`if (!isNumeric(CAR(args)) \|\| !isNumeric(CADR(args)))`
		435	`#endif`
		436	`error("%%*%% requires numeric matrix/vector arguments\n");`
		437
		438	`x = CAR(args);`
		439	`y = CADR(args);`
		440	`xdims = getAttrib(x, R_DimSymbol);`
		441	`ydims = getAttrib(y, R_DimSymbol);`
		442	`ldx = LENGTH(xdims);`
		443	`ldy = LENGTH(ydims);`
		444
		445	`if (ldx != 2 && ldy != 2) {`
		446	`if(PRIMVAL(op) == 0) {`
		447	`nrx = 1;`
		448	`ncx = LENGTH(x);`
		449	`}`
		450	`else {`
		451	`nrx = LENGTH(x);`
		452	`ncx = 1;`
		453	`}`
		454	`nry = LENGTH(y);`
		455	`ncy = 1;`
		456	`}`
		457	`else if (ldx != 2) {`
		458	`nry = INTEGER(ydims)[0];`
		459	`ncy = INTEGER(ydims)[1];`
		460	`nrx = 0;`
		461	`ncx = 0;`
		462	`if(PRIMVAL(op) == 0) {`
		463	`if(LENGTH(x) == nry) {`
		464	`nrx = 1;`
		465	`ncx = LENGTH(x);`
		466	`}`
		467	`}`
		468	`else {`
		469	`if(LENGTH(x) == nry) {`
		470	`nrx = LENGTH(x);`
		471	`ncx = 1;`
		472	`}`
		473	`}`
		474	`}`
		475	`else if (ldy != 2) {`
		476	`nrx = INTEGER(xdims)[0];`
		477	`ncx = INTEGER(xdims)[1];`
		478	`nry = 0;`
		479	`ncy = 0;`
		480	`if(PRIMVAL(op) == 0) {`
		481	`if (LENGTH(y) == ncx) {`
		482	`nry = LENGTH(y);`
		483	`ncy = 1;`
		484	`}`
		485	`}`
		486	`else {`
		487	`if (LENGTH(y) == nrx) {`
		488	`nry = LENGTH(y);`
		489	`ncy = 1;`
		490	`}`
		491	`}`
		492	`}`
		493	`else {`
		494	`nrx = INTEGER(xdims)[0];`
		495	`ncx = INTEGER(xdims)[1];`
		496	`nry = INTEGER(ydims)[0];`
		497	`ncy = INTEGER(ydims)[1];`
		498	`}`
		499
		500	`if(PRIMVAL(op) == 0) {`
		501	`if(ncx != nry)`
		502	`errorcall(call, "non-conformable arguments\n");`
		503	`}`
		504	`else {`
		505	`if(nrx != nry)`
		506	`errorcall(call, "non-conformable arguments\n");`
		507	`}`
		508
		509	`#ifdef COMPLEX_DATA`
		510	`if(isComplex(CAR(args)) \|\| isComplex(CADR(args)))`
		511	`mode = CPLXSXP;`
		512	`else`
		513	`#endif`
		514	`mode = REALSXP;`
		515	`CAR(args) = coerceVector(CAR(args), mode);`
		516	`CADR(args) = coerceVector(CADR(args), mode);`
		517
		518	`if(PRIMVAL(op) == 0) {`
		519	`PROTECT(ans = allocMatrix(mode, nrx, ncy));`
		520	`#ifdef COMPLEX_DATA`
		521	`if(mode == CPLXSXP)`
		522	`cmatprod(COMPLEX(CAR(args)), nrx, ncx, COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));`
		523	`else`
		524	`#endif`
		525	`matprod(REAL(CAR(args)), nrx, ncx, REAL(CADR(args)), nry, ncy, REAL(ans));`
		526	`PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));`
		527	`PROTECT(ydims = getAttrib(CADR(args), R_DimNamesSymbol));`
		528	`if (xdims != R_NilValue \|\| ydims != R_NilValue) {`
		529	`setAttrib(ans, R_DimNamesSymbol, list2(CAR(xdims), CADR(ydims)));`
		530	`}`
		531	`}`
		532	`else {`
		533	`PROTECT(ans = allocMatrix(mode, ncx, ncy));`
		534	`#ifdef COMPLEX_DATA`
		535	`if(mode == CPLXSXP)`
		536	`ccrossprod(COMPLEX(CAR(args)), nrx, ncx, COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));`
		537	`else`
		538	`#endif`
		539	`crossprod(REAL(CAR(args)), nrx, ncx, REAL(CADR(args)), nry, ncy, REAL(ans));`
		540	`PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));`
		541	`PROTECT(ydims = getAttrib(CADR(args), R_DimNamesSymbol));`
		542	`if (xdims != R_NilValue \|\| ydims != R_NilValue) {`
		543	`setAttrib(ans, R_DimNamesSymbol, list2(CADR(xdims), CADR(ydims)));`
		544	`}`
		545	`}`
		546	`UNPROTECT(3);`
		547	`return ans;`
		548	`}`
		549
		550	`SEXP do_transpose(SEXP call, SEXP op, SEXP args, SEXP rho)`
		551	`{`
		552	`SEXP a, r, dims, dn;`
		553	`int i, len, ncol, nrow;`
		554
		555	`checkArity(op, args);`
		556	`a = CAR(args);`
		557
		558	`if(isVector(a)) {`
		559	`dims = getAttrib(a, R_DimSymbol);`
		560	`switch(length(dims)) {`
		561	`case 0:`
		562	`case 1:`
		563	`nrow = len = length(a);`
		564	`ncol = 1;`
		565	`break;`
		566	`case 2:`
		567	`ncol = ncols(a);`
		568	`nrow = nrows(a);`
		569	`len = length(a);`
		570	`break;`
		571	`default:`
		572	`goto not_matrix;`
		573	`}`
		574	`}`
		575	`else if(isList(a)) {`
		576	`dims = getAttrib(a, R_DimSymbol);`
		577	`if(length(dims) == 2) {`
		578	`errorcall(call, "can't transpose list matrices (yet)\n");`
		579	`}`
		580	`else goto not_matrix;`
		581	`}`
		582	`else goto not_matrix;`
		583
		584	`PROTECT(r = allocVector(TYPEOF(a), len));`
		585
		586	`switch (TYPEOF(a)) {`
		587	`case LGLSXP:`
		588	`case FACTSXP:`
		589	`case ORDSXP:`
		590	`case INTSXP:`
		591	`for (i = 0; i < len; i++)`
		592	`INTEGER(r)[i] = INTEGER(a)[(i / ncol) + (i % ncol) * nrow];`
		593	`break;`
		594	`case REALSXP:`
		595	`for (i = 0; i < len; i++)`
		596	`REAL(r)[i] = REAL(a)[(i / ncol) + (i % ncol) * nrow];`
		597	`break;`
		598	`#ifdef COMPLEX_DATA`
		599	`case CPLXSXP:`
		600	`for (i = 0; i < len; i++)`
		601	`COMPLEX(r)[i] = COMPLEX(a)[(i / ncol) + (i % ncol) * nrow];`
		602	`break;`
		603	`#endif`
		604	`case STRSXP:`
		605	`for (i = 0; i < len; i++)`
		606	`STRING(r)[i] = STRING(a)[(i / ncol) + (i % ncol) * nrow];`
		607	`break;`
		608	`}`
		609	`dims = allocVector(INTSXP, 2);`
		610	`INTEGER(dims)[0] = ncol;`
		611	`INTEGER(dims)[1] = nrow;`
		612	`setAttrib(r, R_DimSymbol, dims);`
		613
		614	`if(!isNull(dn = getAttrib(a, R_DimNamesSymbol))) {`
		615	`PROTECT(dn = duplicate(dn));`
		616	`switch(length(dn)) {`
		617	`case 1:`
		618	`PROTECT(dims = allocList(2));`
		619	`CADR(dims) = CAR(dn);`
		620	`setAttrib(r, R_DimNamesSymbol, dims);`
		621	`UNPROTECT(1);`
		622	`break;`
		623	`case 2:`
		624	`dims = CAR(dn);`
		625	`CAR(dn) = CADR(dn);`
		626	`CADR(dn) = dims;`
		627	`setAttrib(r, R_DimNamesSymbol, dn);`
		628	`break;`
		629	`}`
		630	`UNPROTECT(1);`
		631	`}`
		632
		633	`UNPROTECT(1);`
		634	`return r;`
		635
		636	`not_matrix:`
		637	`errorcall(call, "argument is not a matrix\n");`
		638	`}`
		639
		640
		641	`/* swap works by finding for a index i, the position */`
		642	`/* in the array with dimensions dims1 in terms of */`
		643	`/* (i, j, k, l, m, ...); i.e. component-wise position, */`
		644	`/* then permute these to the order of the array with */`
		645	`/* dimensions dims2 and work backwards to an integer */`
		646	`/* offset in this array */`
		647
		648	`static int swap(int ival, SEXP dims1, SEXP dims2, SEXP perm, SEXP ind1, SEXP ind2)`
		649	`{`
		650	`int len, t1, i;`
		651
		652	`len = length(dims1);`
		653	`t1 = ival;`
		654
		655	`for (i = 0; i < len; i++) {`
		656	`INTEGER(ind1)[i] = t1 % INTEGER(dims1)[i];`
		657	`t1 = t1 / INTEGER(dims1)[i];`
		658	`}`
		659
		660	`for (i = 0; i < len; i++)`
		661	`INTEGER(ind2)[i] = INTEGER(ind1)[(INTEGER(perm)[i] - 1)];`
		662
		663	`t1 = INTEGER(ind2)[(len - 1)];`
		664	`for (i = (len - 2); i >= 0; i--) {`
		665	`t1 *= INTEGER(dims2)[i];`
		666	`t1 += INTEGER(ind2)[i];`
		667	`}`
		668	`return t1;`
		669	`}`
		670
		671	`SEXP do_aperm(SEXP call, SEXP op, SEXP args, SEXP rho)`
		672	`{`
		673	`SEXP a, perm, r, dimsa, dimsr, ind1, ind2;`
		674	`int i, j, len;`
		675
		676	`checkArity(op, args);`
		677
		678	`a = CAR(args);`
		679	`PROTECT(dimsa = getAttrib(a, R_DimSymbol));`
		680	`if (dimsa == R_NilValue)`
		681	`error("aperm: invalid first argument, must be an array\n");`
		682
		683	`PROTECT(perm = coerceVector(CADR(args), INTSXP));`
		684	`if (!isVector(perm) \|\| (length(perm) != length(dimsa)))`
		685	`error("aperm: invalid second argument, must be a vector\n");`
		686
		687	`len = length(a);`
		688
		689	`PROTECT(dimsr = allocVector(INTSXP, length(dimsa)));`
		690	`for (i = 0; i < length(dimsa); i++)`
		691	`INTEGER(dimsr)[i] = INTEGER(dimsa)[(INTEGER(perm)[i] - 1)];`
		692
		693	`PROTECT(r = allocVector(TYPEOF(a), len));`
		694	`PROTECT(ind1 = allocVector(INTSXP, LENGTH(dimsa)));`
		695	`PROTECT(ind2 = allocVector(INTSXP, LENGTH(dimsa)));`
		696
		697	`switch (TYPEOF(a)) {`
		698	`case INTSXP:`
		699	`case FACTSXP:`
		700	`case ORDSXP:`
		701	`case LGLSXP:`
		702	`for (i = 0; i < len; i++) {`
		703	`j = swap(i, dimsa, dimsr, perm, ind1, ind2);`
		704	`INTEGER(r)[j] = INTEGER(a)[i];`
		705	`}`
		706	`break;`
		707	`case REALSXP:`
		708	`for (i = 0; i < len; i++) {`
		709	`j = swap(i, dimsa, dimsr, perm, ind1, ind2);`
		710	`REAL(r)[j] = REAL(a)[i];`
		711	`}`
		712	`break;`
		713	`#ifdef COMPLEX_DATA`
		714	`case CPLXSXP:`
		715	`for (i = 0; i < len; i++) {`
		716	`j = swap(i, dimsa, dimsr, perm, ind1, ind2);`
		717	`COMPLEX(r)[j] = COMPLEX(a)[i];`
		718	`}`
		719	`break;`
		720	`#endif`
		721	`case STRSXP:`
		722	`for (i = 0; i < len; i++) {`
		723	`j = swap(i, dimsa, dimsr, perm, ind1, ind2);`
		724	`STRING(r)[j] = STRING(a)[i];`
		725	`}`
		726	`break;`
		727	`default:`
		728	`errorcall(call, "invalid argument\n");`
		729	`}`
		730
		731	`if (INTEGER(CAR(CDDR(args)))[0])`
		732	`setAttrib(r, R_DimSymbol, dimsr);`
		733	`else`
		734	`setAttrib(r, R_DimSymbol, dimsa);`
		735	`UNPROTECT(6);`
		736	`return r;`
		737	`}`

The R Project SVN R

(root)/trunk/src/main/array.c @ 28362 – Rev 2