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

Rev	Author	Line No.	Line
2	r	1	`/*`
2506	maechler	2	`* R : A Computer Language for Statistical Data Analysis`
76535	ripley	3	`* Copyright (C) 1998-2019 The R Core Team`
67965	maechler	4	`* Copyright (C) 2002-2015 The R Foundation`
73160	maechler	5	`* Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka`
2	r	6	`*`
		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`
		9	`* the Free Software Foundation; either version 2 of the License, or`
		10	`* (at your option) any later version.`
		11	`*`
		12	`* This program is distributed in the hope that it will be useful,`
		13	`* but WITHOUT ANY WARRANTY; without even the implied warranty of`
		14	`* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the`
		15	`* GNU General Public License for more details.`
		16	`*`
42307	ripley	17	`* You should have received a copy of the GNU General Public License`
		18	`* along with this program; if not, a copy is available at`
68947	ripley	19	`* https://www.R-project.org/Licenses/`
2	r	20	`*/`
		21
5187	hornik	22	`#ifdef HAVE_CONFIG_H`
7701	hornik	23	`#include <config.h>`
5187	hornik	24	`#endif`
		25
11499	ripley	26	`#include <Defn.h>`
60667	ripley	27	`#include <Internal.h>`
11499	ripley	28	`#include <Rmath.h>`
76540	ripley	29	`#include <R_ext/RS.h> /* for Calloc/Free, F77_CALL */`
		30	`#include <R_ext/BLAS.h>`
68663	luke	31	`#include <R_ext/Itermacros.h>`
1839	ihaka	32
68663	luke	33	`#include "duplicate.h"`
		34
72294	kalibera	35	`#include <complex.h>`
		36	`#include "Rcomplex.h" /* toC99 */`
		37
26467	maechler	38	`/* "GetRowNames" and "GetColNames" are utility routines which`
		39	`* locate and return the row names and column names from the`
		40	`* dimnames attribute of a matrix. They are useful because`
		41	`* old versions of R used pair-based lists for dimnames`
		42	`* whereas recent versions use vector based lists.`
1881	ihaka	43
26467	maechler	44	`* These are now very old, plus`
		45	* ``When the "dimnames" attribute is
		46	`* grabbed off an array it is always adjusted to be a vector.''`
53676	ripley	47
		48	`They are used in bind.c and subset.c, and advertised in Rinternals.h`
26467	maechler	49	`*/`
1839	ihaka	50	`SEXP GetRowNames(SEXP dimnames)`
		51	`{`
		52	`if (TYPEOF(dimnames) == VECSXP)`
10172	luke	53	`return VECTOR_ELT(dimnames, 0);`
1839	ihaka	54	`else`
1858	ihaka	55	`return R_NilValue;`
2213	maechler	56	`}`
1839	ihaka	57
		58	`SEXP GetColNames(SEXP dimnames)`
		59	`{`
		60	`if (TYPEOF(dimnames) == VECSXP)`
10172	luke	61	`return VECTOR_ELT(dimnames, 1);`
1839	ihaka	62	`else`
2213	maechler	63	`return R_NilValue;`
1839	ihaka	64	`}`
		65
36990	ripley	66	`SEXP attribute_hidden do_matrix(SEXP call, SEXP op, SEXP args, SEXP rho)`
2	r	67	`{`
53676	ripley	68	`SEXP vals, ans, snr, snc, dimnames;`
59052	ripley	69	`int nr = 1, nc = 1, byrow, miss_nr, miss_nc;`
		70	`R_xlen_t lendat;`
2	r	71
1820	ihaka	72	`checkArity(op, args);`
53697	ripley	73	`vals = CAR(args); args = CDR(args);`
		74	`switch(TYPEOF(vals)) {`
		75	`case LGLSXP:`
		76	`case INTSXP:`
		77	`case REALSXP:`
		78	`case CPLXSXP:`
		79	`case STRSXP:`
		80	`case RAWSXP:`
		81	`case EXPRSXP:`
		82	`case VECSXP:`
		83	`break;`
		84	`default:`
61566	ripley	85	`error(_("'data' must be of a vector type, was '%s'"),`
		86	`type2char(TYPEOF(vals)));`
53697	ripley	87	`}`
59052	ripley	88	`lendat = XLENGTH(vals);`
53697	ripley	89	`snr = CAR(args); args = CDR(args);`
		90	`snc = CAR(args); args = CDR(args);`
		91	`byrow = asLogical(CAR(args)); args = CDR(args);`
31894	ripley	92	`if (byrow == NA_INTEGER)`
44512	ripley	93	`error(_("invalid '%s' argument"), "byrow");`
53697	ripley	94	`dimnames = CAR(args);`
53900	ripley	95	`args = CDR(args);`
		96	`miss_nr = asLogical(CAR(args)); args = CDR(args);`
		97	`miss_nc = asLogical(CAR(args));`
2	r	98
53697	ripley	99	`if (!miss_nr) {`
		100	`if (!isNumeric(snr)) error(_("non-numeric matrix extent"));`
		101	`nr = asInteger(snr);`
		102	`if (nr == NA_INTEGER)`
		103	`error(_("invalid 'nrow' value (too large or NA)"));`
		104	`if (nr < 0)`
		105	`error(_("invalid 'nrow' value (< 0)"));`
		106	`}`
		107	`if (!miss_nc) {`
		108	`if (!isNumeric(snc)) error(_("non-numeric matrix extent"));`
		109	`nc = asInteger(snc);`
		110	`if (nc == NA_INTEGER)`
		111	`error(_("invalid 'ncol' value (too large or NA)"));`
		112	`if (nc < 0)`
		113	`error(_("invalid 'ncol' value (< 0)"));`
		114	`}`
59090	ripley	115	`if (miss_nr && miss_nc) {`
		116	`if (lendat > INT_MAX) error("data is too long");`
		117	`nr = (int) lendat;`
		118	`} else if (miss_nr) {`
		119	`if (lendat > (double) nc * INT_MAX) error("data is too long");`
64246	ripley	120	`// avoid division by zero`
		121	`if (nc == 0) {`
		122	`if (lendat) error(_("nc = 0 for non-null data"));`
		123	`else nr = 0;`
		124	`} else`
		125	`nr = (int) ceil((double) lendat / (double) nc);`
59090	ripley	126	`} else if (miss_nc) {`
		127	`if (lendat > (double) nr * INT_MAX) error("data is too long");`
64246	ripley	128	`// avoid division by zero`
		129	`if (nr == 0) {`
		130	`if (lendat) error(_("nr = 0 for non-null data"));`
		131	`else nc = 0;`
		132	`} else`
		133	`nc = (int) ceil((double) lendat / (double) nr);`
59090	ripley	134	`}`
2	r	135
59090	ripley	136	`if(lendat > 0) {`
59052	ripley	137	`R_xlen_t nrc = (R_xlen_t) nr * nc;`
		138	`if (lendat > 1 && nrc % lendat != 0) {`
28014	ripley	139	`if (((lendat > nr) && (lendat / nr) * nr != lendat) \|\|`
		140	`((lendat < nr) && (nr / lendat) * lendat != nr))`
41686	ripley	141	`warning(_("data length [%d] is not a sub-multiple or multiple of the number of rows [%d]"), lendat, nr);`
28014	ripley	142	`else if (((lendat > nc) && (lendat / nc) * nc != lendat) \|\|`
		143	`((lendat < nc) && (nc / lendat) * lendat != nc))`
41686	ripley	144	`warning(_("data length [%d] is not a sub-multiple or multiple of the number of columns [%d]"), lendat, nc);`
28014	ripley	145	`}`
59052	ripley	146	`else if ((lendat > 1) && (nrc == 0)){`
32860	ripley	147	`warning(_("data length exceeds size of matrix"));`
2	r	148	`}`
28014	ripley	149	`}`
2	r	150
61696	maechler	151	`#ifndef LONG_VECTOR_SUPPORT`
		152	`if ((double)nr * (double)nc > INT_MAX)`
41680	ripley	153	`error(_("too many elements specified"));`
59052	ripley	154	`#endif`
28349	ripley	155
43007	ripley	156	`PROTECT(ans = allocMatrix(TYPEOF(vals), nr, nc));`
28014	ripley	157	`if(lendat) {`
		158	`if (isVector(vals))`
43007	ripley	159	`copyMatrix(ans, vals, byrow);`
28014	ripley	160	`else`
43007	ripley	161	`copyListMatrix(ans, vals, byrow);`
28014	ripley	162	`} else if (isVector(vals)) { /* fill with NAs */`
60540	urbaneks	163	`R_xlen_t N = (R_xlen_t) nr * nc, i;`
28014	ripley	164	`switch(TYPEOF(vals)) {`
		165	`case STRSXP:`
60540	urbaneks	166	`for (i = 0; i < N; i++)`
		167	`SET_STRING_ELT(ans, i, NA_STRING);`
28014	ripley	168	`break;`
		169	`case LGLSXP:`
60540	urbaneks	170	`for (i = 0; i < N; i++)`
		171	`LOGICAL(ans)[i] = NA_LOGICAL;`
28014	ripley	172	`break;`
		173	`case INTSXP:`
60540	urbaneks	174	`for (i = 0; i < N; i++)`
		175	`INTEGER(ans)[i] = NA_INTEGER;`
28014	ripley	176	`break;`
		177	`case REALSXP:`
60540	urbaneks	178	`for (i = 0; i < N; i++)`
		179	`REAL(ans)[i] = NA_REAL;`
28014	ripley	180	`break;`
		181	`case CPLXSXP:`
		182	`{`
		183	`Rcomplex na_cmplx;`
		184	`na_cmplx.r = NA_REAL;`
		185	`na_cmplx.i = 0;`
60540	urbaneks	186	`for (i = 0; i < N; i++)`
		187	`COMPLEX(ans)[i] = na_cmplx;`
28014	ripley	188	`}`
		189	`break;`
29902	ripley	190	`case RAWSXP:`
69111	ripley	191	`if (N) memset(RAW(ans), 0, N);`
29902	ripley	192	`break;`
32350	maechler	193	`default:`
31875	ripley	194	`/* don't fill with anything */`
		195	`;`
28014	ripley	196	`}`
		197	`}`
69786	maechler	198	`if(!isNull(dimnames) && length(dimnames) > 0)`
47049	ripley	199	`ans = dimnamesgets(ans, dimnames);`
1820	ihaka	200	`UNPROTECT(1);`
43007	ripley	201	`return ans;`
2	r	202	`}`
		203
		204
		205	`SEXP allocMatrix(SEXPTYPE mode, int nrow, int ncol)`
		206	`{`
1820	ihaka	207	`SEXP s, t;`
59052	ripley	208	`R_xlen_t n;`
2	r	209
1820	ihaka	210	`if (nrow < 0 \|\| ncol < 0)`
32860	ripley	211	`error(_("negative extents to matrix"));`
59052	ripley	212	`#ifndef LONG_VECTOR_SUPPORT`
28349	ripley	213	`if ((double)nrow * (double)ncol > INT_MAX)`
32860	ripley	214	`error(_("allocMatrix: too many elements specified"));`
59052	ripley	215	`#endif`
		216	`n = ((R_xlen_t) nrow) * ncol;`
1820	ihaka	217	`PROTECT(s = allocVector(mode, n));`
		218	`PROTECT(t = allocVector(INTSXP, 2));`
		219	`INTEGER(t)[0] = nrow;`
		220	`INTEGER(t)[1] = ncol;`
		221	`setAttrib(s, R_DimSymbol, t);`
		222	`UNPROTECT(2);`
		223	`return s;`
2	r	224	`}`
		225
42325	bates	226	`/**`
		227	`* Allocate a 3-dimensional array`
		228	`*`
		229	`* @param mode The R mode (e.g. INTSXP)`
		230	`* @param nrow number of rows`
		231	`* @param ncol number of columns`
		232	`* @param nface number of faces`
		233	`*`
		234	`* @return A 3-dimensional array of the indicated dimensions and mode`
		235	`*/`
		236	`SEXP alloc3DArray(SEXPTYPE mode, int nrow, int ncol, int nface)`
		237	`{`
		238	`SEXP s, t;`
59052	ripley	239	`R_xlen_t n;`
1130	maechler	240
42325	bates	241	`if (nrow < 0 \|\| ncol < 0 \|\| nface < 0)`
		242	`error(_("negative extents to 3D array"));`
59052	ripley	243	`#ifndef LONG_VECTOR_SUPPORT`
42325	bates	244	`if ((double)nrow * (double)ncol * (double)nface > INT_MAX)`
75324	maechler	245	`error(_("'alloc3DArray': too many elements specified"));`
59052	ripley	246	`#endif`
		247	`n = ((R_xlen_t) nrow) * ncol * nface;`
42325	bates	248	`PROTECT(s = allocVector(mode, n));`
		249	`PROTECT(t = allocVector(INTSXP, 3));`
		250	`INTEGER(t)[0] = nrow;`
		251	`INTEGER(t)[1] = ncol;`
		252	`INTEGER(t)[2] = nface;`
		253	`setAttrib(s, R_DimSymbol, t);`
		254	`UNPROTECT(2);`
		255	`return s;`
		256	`}`
		257
		258
2	r	259	`SEXP allocArray(SEXPTYPE mode, SEXP dims)`
		260	`{`
1820	ihaka	261	`SEXP array;`
59052	ripley	262	`int i;`
59090	ripley	263	`R_xlen_t n = 1;`
75325	ripley	264	`#ifndef LONG_VECTOR_SUPPORT`
59090	ripley	265	`double dn = 1;`
75325	ripley	266	`#endif`
2	r	267
28349	ripley	268	`for (i = 0; i < LENGTH(dims); i++) {`
75324	maechler	269	`#ifndef LONG_VECTOR_SUPPORT`
28349	ripley	270	`dn *= INTEGER(dims)[i];`
		271	`if(dn > INT_MAX)`
60844	ripley	272	`error(_("'allocArray': too many elements specified by 'dims'"));`
59052	ripley	273	`#endif`
28349	ripley	274	`n *= INTEGER(dims)[i];`
		275	`}`
2	r	276
1820	ihaka	277	`PROTECT(dims = duplicate(dims));`
		278	`PROTECT(array = allocVector(mode, n));`
		279	`setAttrib(array, R_DimSymbol, dims);`
		280	`UNPROTECT(2);`
		281	`return array;`
2	r	282	`}`
		283
1820	ihaka	284	`/* DropDims strips away redundant dimensioning information. */`
		285	`/* If there is an appropriate dimnames attribute the correct */`
		286	`/* element is extracted and attached to the vector as a names */`
		287	`/* attribute. Note that this function mutates x. */`
		288	`/* Duplication should occur before this is called. */`
2	r	289
		290	`SEXP DropDims(SEXP x)`
		291	`{`
43878	ripley	292	`SEXP dims, dimnames, newnames = R_NilValue;`
1839	ihaka	293	`int i, n, ndims;`
2	r	294
1820	ihaka	295	`PROTECT(x);`
		296	`dims = getAttrib(x, R_DimSymbol);`
2	r	297
1820	ihaka	298	`/* Check that dropping will actually do something. */`
		299	`/* (1) Check that there is a "dim" attribute. */`
2	r	300
1839	ihaka	301	`if (dims == R_NilValue) {`
68092	luke	302	`UNPROTECT(1); /* x */`
1820	ihaka	303	`return x;`
		304	`}`
1839	ihaka	305	`ndims = LENGTH(dims);`
69813	maechler	306	`int *dim = INTEGER(dims); // used several times`
2	r	307
1839	ihaka	308	`/* (2) Check whether there are redundant extents */`
1820	ihaka	309	`n = 0;`
1839	ihaka	310	`for (i = 0; i < ndims; i++)`
69813	maechler	311	`if (dim[i] != 1) n++;`
1839	ihaka	312	`if (n == ndims) {`
68092	luke	313	`UNPROTECT(1); /* x */`
1820	ihaka	314	`return x;`
		315	`}`
2	r	316
68092	luke	317	`PROTECT(dimnames = getAttrib(x, R_DimNamesSymbol));`
1839	ihaka	318	`if (n <= 1) {`
43878	ripley	319	`/* We have reduced to a vector result.`
		320	`If that has length one, it is ambiguous which dimnames to use,`
		321	`so use it if there is only one (as from R 2.7.0).`
		322	`*/`
1820	ihaka	323	`if (dimnames != R_NilValue) {`
59052	ripley	324	`if(XLENGTH(x) != 1) {`
43878	ripley	325	`for (i = 0; i < LENGTH(dims); i++) {`
69813	maechler	326	`if (dim[i] != 1) {`
10172	luke	327	`newnames = VECTOR_ELT(dimnames, i);`
1839	ihaka	328	`break;`
		329	`}`
2	r	330	`}`
43878	ripley	331	`} else { /* drop all dims: keep names if unambiguous */`
		332	`int cnt;`
		333	`for(i = 0, cnt = 0; i < LENGTH(dims); i++)`
		334	`if(VECTOR_ELT(dimnames, i) != R_NilValue) cnt++;`
		335	`if(cnt == 1)`
		336	`for (i = 0; i < LENGTH(dims); i++) {`
		337	`newnames = VECTOR_ELT(dimnames, i);`
		338	`if(newnames != R_NilValue) break;`
1839	ihaka	339	`}`
1820	ihaka	340	`}`
2	r	341	`}`
1820	ihaka	342	`PROTECT(newnames);`
		343	`setAttrib(x, R_DimNamesSymbol, R_NilValue);`
		344	`setAttrib(x, R_DimSymbol, R_NilValue);`
		345	`setAttrib(x, R_NamesSymbol, newnames);`
48062	jmc	346	`/* FIXME: the following is desirable, but pointless as long as`
		347	`subset.c & others have a contrary version that leaves the`
		348	`S4 class in, incorrectly, in the case of vectors. JMC`
		349	`3/3/09 */`
68923	ripley	350	`/* if(IS_S4_OBJECT(x)) {/\* no longer valid subclass of array or`
		351	`matrix \/ /`
		352	`/* setAttrib(x, R_ClassSymbol, R_NilValue); */`
		353	`/* UNSET_S4_OBJECT(x); */`
		354	`/* } */`
68092	luke	355	`UNPROTECT(1); /* newnames */`
43878	ripley	356	`} else {`
69813	maechler	357	`// We have a lower dimensional array, and n == length(newdims)`
6994	pd	358	`SEXP newdims, dnn, newnamesnames = R_NilValue;`
68092	luke	359	`PROTECT(dnn = getAttrib(dimnames, R_NamesSymbol));`
1820	ihaka	360	`PROTECT(newdims = allocVector(INTSXP, n));`
1858	ihaka	361	`for (i = 0, n = 0; i < ndims; i++)`
69813	maechler	362	`if (dim[i] != 1)`
		363	`INTEGER(newdims)[n++] = dim[i];`
		364	`if(!isNull(getAttrib(dims, R_NamesSymbol))) {`
		365	`SEXP nms_d = getAttrib(dims, R_NamesSymbol),`
		366	`new_nms = PROTECT(allocVector(STRSXP, n));`
		367	`for (i = 0, n = 0; i < ndims; i++)`
		368	`if (dim[i] != 1)`
		369	`SET_STRING_ELT(new_nms, n++, STRING_ELT(nms_d, i));`
		370	`setAttrib(newdims, R_NamesSymbol, new_nms);`
		371	`UNPROTECT(1);`
		372	`}`
		373	`Rboolean havenames = FALSE;`
6994	pd	374	`if (!isNull(dimnames)) {`
1858	ihaka	375	`for (i = 0; i < ndims; i++)`
69813	maechler	376	`if (dim[i] != 1 &&`
10172	luke	377	`VECTOR_ELT(dimnames, i) != R_NilValue)`
69813	maechler	378	`havenames = TRUE;`
1858	ihaka	379	`if (havenames) {`
2715	pd	380	`PROTECT(newnames = allocVector(VECSXP, n));`
6994	pd	381	`PROTECT(newnamesnames = allocVector(STRSXP, n));`
		382	`for (i = 0, n = 0; i < ndims; i++) {`
69813	maechler	383	`if (dim[i] != 1) {`
6994	pd	384	`if(!isNull(dnn))`
10172	luke	385	`SET_STRING_ELT(newnamesnames, n,`
		386	`STRING_ELT(dnn, i));`
		387	`SET_VECTOR_ELT(newnames, n++, VECTOR_ELT(dimnames, i));`
6994	pd	388	`}`
1858	ihaka	389	`}`
1820	ihaka	390	`}`
1858	ihaka	391	`else dimnames = R_NilValue;`
1839	ihaka	392	`}`
1820	ihaka	393	`setAttrib(x, R_DimNamesSymbol, R_NilValue);`
		394	`setAttrib(x, R_DimSymbol, newdims);`
69813	maechler	395	`if (havenames)`
2715	pd	396	`{`
6994	pd	397	`if(!isNull(dnn))`
		398	`setAttrib(newnames, R_NamesSymbol, newnamesnames);`
2506	maechler	399	`setAttrib(x, R_DimNamesSymbol, newnames);`
68092	luke	400	`UNPROTECT(2); /* newnamesnames, newnames */`
2715	pd	401	`}`
68092	luke	402	`UNPROTECT(2); /* newdims, dnn */`
1820	ihaka	403	`}`
68092	luke	404	`UNPROTECT(2); /* dimnames, x */`
1820	ihaka	405	`return x;`
2	r	406	`}`
		407
36990	ripley	408	`SEXP attribute_hidden do_drop(SEXP call, SEXP op, SEXP args, SEXP rho)`
2	r	409	`{`
1820	ihaka	410	`SEXP x, xdims;`
		411	`int i, n, shorten;`
2	r	412
1820	ihaka	413	`checkArity(op, args);`
		414	`x = CAR(args);`
1839	ihaka	415	`if ((xdims = getAttrib(x, R_DimSymbol)) != R_NilValue) {`
1820	ihaka	416	`n = LENGTH(xdims);`
		417	`shorten = 0;`
1839	ihaka	418	`for (i = 0; i < n; i++)`
		419	`if (INTEGER(xdims)[i] == 1) shorten = 1;`
		420	`if (shorten) {`
75534	luke	421	`if (MAYBE_REFERENCED(x)) x = R_duplicate_attr(x);`
1820	ihaka	422	`x = DropDims(x);`
2	r	423	`}`
1820	ihaka	424	`}`
		425	`return x;`
2	r	426	`}`
		427
1820	ihaka	428	`/* Length of Primitive Objects */`
2	r	429
36990	ripley	430	`SEXP attribute_hidden do_length(SEXP call, SEXP op, SEXP args, SEXP rho)`
2	r	431	`{`
40176	ripley	432	`checkArity(op, args);`
51267	ripley	433	`check1arg(args, call, "x");`
15783	rgentlem	434
60156	ripley	435	`SEXP x = CAR(args), ans;`
59185	ripley	436
60156	ripley	437	`if (isObject(x) &&`
59239	ripley	438	`DispatchOrEval(call, op, "length", args, rho, &ans, 0, 1)) {`
		439	`if (length(ans) == 1 && TYPEOF(ans) == REALSXP) {`
59241	ripley	440	`double d = REAL(ans)[0];`
72383	morgan	441	`if (R_FINITE(d) && d >= 0. && d <= INT_MAX && floor(d) == d) {`
		442	`PROTECT(ans);`
		443	`ans = coerceVector(ans, INTSXP);`
		444	`UNPROTECT(1);`
		445	`return(ans);`
		446	`}`
59239	ripley	447	`}`
59224	ripley	448	`return(ans);`
59239	ripley	449	`}`
17515	ripley	450
67965	maechler	451
59185	ripley	452	`#ifdef LONG_VECTOR_SUPPORT`
60132	ripley	453	`// or use IS_LONG_VEC`
59185	ripley	454	`R_xlen_t len = xlength(x);`
60132	ripley	455	`if (len > INT_MAX) return ScalarReal((double) len);`
		456	`#endif`
59185	ripley	457	`return ScalarInteger(length(x));`
2	r	458	`}`
		459
70308	lawrence	460	`R_len_t attribute_hidden dispatch_length(SEXP x, SEXP call, SEXP rho) {`
		461	`R_xlen_t len = dispatch_xlength(x, call, rho);`
70312	ripley	462	`#ifdef LONG_VECTOR_SUPPORT`
		463	`if (len > INT_MAX) return R_BadLongVector(x, __FILE__, __LINE__);`
		464	`#endif`
		465	`return (R_len_t) len;`
70308	lawrence	466	`}`
		467
		468	`R_xlen_t attribute_hidden dispatch_xlength(SEXP x, SEXP call, SEXP rho) {`
67750	lawrence	469	`static SEXP length_op = NULL;`
69448	lawrence	470	`if (isObject(x)) {`
		471	`SEXP len, args;`
		472	`if (length_op == NULL)`
		473	`length_op = R_Primitive("length");`
		474	`PROTECT(args = list1(x));`
		475	`if (DispatchOrEval(call, length_op, "length", args, rho, &len, 0, 1)) {`
		476	`UNPROTECT(1);`
		477	`return (R_xlen_t)`
		478	`(TYPEOF(len) == REALSXP ? REAL(len)[0] : asInteger(len));`
		479	`}`
		480	`UNPROTECT(1);`
67750	lawrence	481	`}`
69642	maechler	482	`return(xlength(x));`
67750	lawrence	483	`}`
2	r	484
69448	lawrence	485	`// auxiliary for do_lengths_*(), i.e., R's lengths()`
		486	`static R_xlen_t getElementLength(SEXP x, R_xlen_t i, SEXP call, SEXP rho) {`
74030	kalibera	487	`SEXP x_elt;`
		488	`R_xlen_t ans;`
		489
74029	kalibera	490	`PROTECT(x_elt = dispatch_subset2(x, i, call, rho));`
		491	`ans = dispatch_xlength(x_elt, call, rho);`
		492	`UNPROTECT(1); /* x_elt */`
		493	`return ans;`
69448	lawrence	494	`}`
		495
68033	ripley	496	`#ifdef LONG_VECTOR_SUPPORT`
67750	lawrence	497	`static SEXP do_lengths_long(SEXP x, SEXP call, SEXP rho)`
		498	`{`
		499	`SEXP ans;`
		500	`R_xlen_t x_len, i;`
		501	`double *ans_elt;`
67965	maechler	502
70308	lawrence	503	`x_len = dispatch_xlength(x, call, rho);`
67750	lawrence	504	`PROTECT(ans = allocVector(REALSXP, x_len));`
68033	ripley	505	`for (i = 0, ans_elt = REAL(ans); i < x_len; i++, ans_elt++)`
72154	ripley	506	`*ans_elt = (double) getElementLength(x, i, call, rho);`
67750	lawrence	507	`UNPROTECT(1);`
		508	`return ans;`
		509	`}`
68033	ripley	510	`#endif`
67750	lawrence	511
		512	`SEXP attribute_hidden do_lengths(SEXP call, SEXP op, SEXP args, SEXP rho)`
		513	`{`
69326	luke	514	`checkArity(op, args);`
67750	lawrence	515	`SEXP x = CAR(args), ans;`
		516	`R_xlen_t x_len, i;`
		517	`int *ans_elt;`
67965	maechler	518	`int useNames = asLogical(CADR(args));`
		519	`if (useNames == NA_LOGICAL)`
70305	lawrence	520	`error(_("invalid '%s' value"), "use.names");`
70219	lawrence	521
		522	`if (DispatchOrEval(call, op, "lengths", args, rho, &ans, 0, 1))`
		523	`return(ans);`
		524
69448	lawrence	525	`Rboolean isList = isVectorList(x) \|\| isS4(x);`
68292	maechler	526	`if(!isList) switch(TYPEOF(x)) {`
		527	`case NILSXP:`
		528	`case CHARSXP:`
		529	`case LGLSXP:`
		530	`case INTSXP:`
		531	`case REALSXP:`
		532	`case CPLXSXP:`
		533	`case STRSXP:`
		534	`case RAWSXP:`
		535	`break;`
		536	`default:`
		537	`error(_("'%s' must be a list or atomic vector"), "x");`
		538	`}`
70308	lawrence	539	`x_len = dispatch_xlength(x, call, rho);`
67750	lawrence	540	`PROTECT(ans = allocVector(INTSXP, x_len));`
68292	maechler	541	`if(isList) {`
		542	`for (i = 0, ans_elt = INTEGER(ans); i < x_len; i++, ans_elt++) {`
		543	`R_xlen_t x_elt_len = getElementLength(x, i, call, rho);`
67750	lawrence	544	`#ifdef LONG_VECTOR_SUPPORT`
68292	maechler	545	`if (x_elt_len > INT_MAX) {`
		546	`ans = do_lengths_long(x, call, rho);`
		547	`UNPROTECT(1);`
		548	`PROTECT(ans);`
		549	`break;`
		550	`}`
67750	lawrence	551	`#endif`
68292	maechler	552	`*ans_elt = (int)x_elt_len;`
		553	`}`
		554	`} else { // atomic: every element has length 1`
		555	`for (i = 0, ans_elt = INTEGER(ans); i < x_len; i++, ans_elt++)`
		556	`*ans_elt = 1;`
67750	lawrence	557	`}`
71907	lawrence	558	`SEXP dim = getAttrib(x, R_DimSymbol);`
		559	`if(!isNull(dim)) {`
		560	`setAttrib(ans, R_DimSymbol, dim);`
		561	`}`
67965	maechler	562	`if(useNames) {`
		563	`SEXP names = getAttrib(x, R_NamesSymbol);`
		564	`if(!isNull(names)) setAttrib(ans, R_NamesSymbol, names);`
71907	lawrence	565	`SEXP dimnames = getAttrib(x, R_DimNamesSymbol);`
71927	lawrence	566	`if(!isNull(dimnames)) setAttrib(ans, R_DimNamesSymbol, dimnames);`
67965	maechler	567	`}`
68187	luke	568	`UNPROTECT(1);`
67750	lawrence	569	`return ans;`
		570	`}`
		571
36990	ripley	572	`SEXP attribute_hidden do_rowscols(SEXP call, SEXP op, SEXP args, SEXP rho)`
2	r	573	`{`
41712	ripley	574	`checkArity(op, args);`
75324	maechler	575	`SEXP dim = CAR(args);`
		576	`int nprot = 0;`
		577	`if (!isInteger(dim)) {`
		578	`PROTECT(dim = coerceVector(dim, INTSXP)); nprot++;`
		579	`}`
		580	`if (LENGTH(dim) != 2)`
60844	ripley	581	`error(_("a matrix-like object is required as argument to '%s'"),`
		582	`(PRIMVAL(op) == 2) ? "col" : "row");`
2	r	583
75324	maechler	584	`int nr = INTEGER(dim)[0],`
		585	`nc = INTEGER(dim)[1];`
		586	`if(nprot) UNPROTECT(nprot);`
2	r	587
75324	maechler	588	`SEXP ans = allocMatrix(INTSXP, nr, nc);`
2	r	589
59052	ripley	590	`R_xlen_t NR = nr;`
1820	ihaka	591	`switch (PRIMVAL(op)) {`
75324	maechler	592	`case 1: // row() & .row()`
		593	`for (int i = 0; i < nr; i++)`
		594	`for (int j = 0; j < nc; j++)`
59052	ripley	595	`INTEGER(ans)[i + j * NR] = i + 1;`
1820	ihaka	596	`break;`
75324	maechler	597	`case 2: // col() & .col()`
		598	`for (int i = 0; i < nr; i++)`
		599	`for (int j = 0; j < nc; j++)`
59052	ripley	600	`INTEGER(ans)[i + j * NR] = j + 1;`
1820	ihaka	601	`break;`
		602	`}`
		603	`return ans;`
2	r	604	`}`
		605
72020	kalibera	606	`/*`
		607	`Whenever vector x contains NaN or Inf (or -Inf), the function returns TRUE.`
		608	`It can be imprecise: it can return TRUE in other cases as well.`
		609
		610	`A precise version of the function could be implemented as`
		611
		612	`for (R_xlen_t i = 0; i < n; i++)`
		613	`if (!R_FINITE(x[i])) return TRUE;`
		614	`return FALSE;`
		615
		616	`The present version is imprecise, but faster.`
		617	`*/`
		618	`static Rboolean mayHaveNaNOrInf(double *x, R_xlen_t n)`
		619	`{`
		620	`if ((n&1) != 0 && !R_FINITE(x[0]))`
		621	`return TRUE;`
72059	kalibera	622	`for (R_xlen_t i = n&1; i < n; i += 2)`
72020	kalibera	623	`/* A precise version could use this condition:`
		624	`*`
		625	`* !R_FINITE(x[i]+x[i+1]) && (!R_FINITE(x[i]) \|\| !R_FINITE(x[i+1]))`
		626	`*`
		627	`* The present imprecise version has been found to be faster`
		628	`* with GCC and ICC in the common case when the sum of the two`
		629	`* values is always finite.`
		630	`*`
		631	`* The present version is imprecise because the sum of two very`
		632	`* large finite values (e.g. 1e308) may be infinite.`
		633	`*/`
		634	`if (!R_FINITE(x[i]+x[i+1]))`
		635	`return TRUE;`
		636	`return FALSE;`
		637	`}`
		638
72142	kalibera	639	`/*`
		640	`This is an experimental version that has been observed to run fast on some`
		641	`SIMD hardware with GCC and ICC.`
		642	`Note that the OpenMP reduction assumes associativity of addition, which is`
		643	`safe here, because the result is only used for an imprecise test for`
		644	`the presence of NaN and Inf values.`
		645	`*/`
		646	`static Rboolean mayHaveNaNOrInf_simd(double *x, R_xlen_t n)`
		647	`{`
		648	`double s = 0;`
72179	kalibera	649	`/* SIMD reduction is supported since OpenMP 4.0. The value of _OPENMP is`
		650	`unreliable in some compilers, so we depend on HAVE_OPENMP_SIMDRED,`
		651	`which is normally set by configure based on a test. */`
		652	`/* _OPENMP >= 201307 */`
		653	`#if defined(_OPENMP) && HAVE_OPENMP_SIMDRED`
72142	kalibera	654	`#pragma omp simd reduction(+:s)`
72174	kalibera	655	`#endif`
72142	kalibera	656	`for (R_xlen_t i = 0; i < n; i++)`
		657	`s += x[i];`
		658	`return !R_FINITE(s);`
		659	`}`
		660
72083	kalibera	661	`static Rboolean cmayHaveNaNOrInf(Rcomplex *x, R_xlen_t n)`
		662	`{`
		663	`/* With HAVE_FORTRAN_DOUBLE_COMPLEX set, it should be clear that`
		664	`Rcomplex has no padding, so we could probably use mayHaveNaNOrInf,`
		665	`but better safe than sorry... */`
		666	`if ((n&1) != 0 && (!R_FINITE(x[0].r) \|\| !R_FINITE(x[0].i)))`
		667	`return TRUE;`
		668	`for (R_xlen_t i = n&1; i < n; i += 2)`
		669	`if (!R_FINITE(x[i].r+x[i].i+x[i+1].r+x[i+1].i))`
		670	`return TRUE;`
		671	`return FALSE;`
		672	`}`
		673
72142	kalibera	674	`/* experimental version for SIMD hardware (see also mayHaveNaNOrInf_simd) */`
		675	`static Rboolean cmayHaveNaNOrInf_simd(Rcomplex *x, R_xlen_t n)`
72059	kalibera	676	`{`
72142	kalibera	677	`double s = 0;`
72179	kalibera	678	`/* _OPENMP >= 201307 - see mayHaveNaNOrInf_simd */`
		679	`#if defined(_OPENMP) && HAVE_OPENMP_SIMDRED`
72142	kalibera	680	`#pragma omp simd reduction(+:s)`
72174	kalibera	681	`#endif`
72142	kalibera	682	`for (R_xlen_t i = 0; i < n; i++) {`
		683	`s += x[i].r;`
		684	`s += x[i].i;`
		685	`}`
		686	`return !R_FINITE(s);`
		687	`}`
		688
		689	`static void internal_matprod(double *x, int nrx, int ncx,`
		690	`double y, int nry, int ncy, double z)`
		691	`{`
72059	kalibera	692	`LDOUBLE sum;`
72142	kalibera	693	`#define MATPROD_BODY \`
		694	`R_xlen_t NRX = nrx, NRY = nry; \`
		695	`for (int i = 0; i < nrx; i++) \`
		696	`for (int k = 0; k < ncy; k++) { \`
		697	`sum = 0.0; \`
		698	`for (int j = 0; j < ncx; j++) \`
		699	`sum += x[i + j * NRX] * y[j + k * NRY]; \`
72153	ripley	700	`z[i + k * NRX] = (double) sum; \`
72059	kalibera	701	`}`
72142	kalibera	702	`MATPROD_BODY;`
72059	kalibera	703	`}`
		704
72142	kalibera	705	`static void simple_matprod(double *x, int nrx, int ncx,`
		706	`double y, int nry, int ncy, double z)`
72059	kalibera	707	`{`
72142	kalibera	708	`double sum;`
		709	`MATPROD_BODY;`
		710	`}`
		711
		712	`static void internal_crossprod(double *x, int nrx, int ncx,`
		713	`double y, int nry, int ncy, double z)`
		714	`{`
72059	kalibera	715	`LDOUBLE sum;`
72142	kalibera	716	`#define CROSSPROD_BODY \`
		717	`R_xlen_t NRX = nrx, NRY = nry, NCX = ncx; \`
		718	`for (int i = 0; i < ncx; i++) \`
		719	`for (int k = 0; k < ncy; k++) { \`
		720	`sum = 0.0; \`
		721	`for (int j = 0; j < nrx; j++) \`
		722	`sum += x[j + i * NRX] * y[j + k * NRY]; \`
72153	ripley	723	`z[i + k * NCX] = (double) sum; \`
72059	kalibera	724	`}`
72142	kalibera	725	`CROSSPROD_BODY;`
72059	kalibera	726	`}`
		727
72142	kalibera	728	`static void simple_crossprod(double *x, int nrx, int ncx,`
		729	`double y, int nry, int ncy, double z)`
72059	kalibera	730	`{`
72142	kalibera	731	`double sum;`
		732	`CROSSPROD_BODY;`
		733	`}`
		734
		735	`static void internal_tcrossprod(double *x, int nrx, int ncx,`
		736	`double y, int nry, int ncy, double z)`
		737	`{`
72059	kalibera	738	`LDOUBLE sum;`
72142	kalibera	739	`#define TCROSSPROD_BODY \`
		740	`R_xlen_t NRX = nrx, NRY = nry; \`
		741	`for (int i = 0; i < nrx; i++) \`
		742	`for (int k = 0; k < nry; k++) { \`
		743	`sum = 0.0; \`
		744	`for (int j = 0; j < ncx; j++) \`
		745	`sum += x[i + j * NRX] * y[k + j * NRY]; \`
72153	ripley	746	`z[i + k * NRX] = (double) sum; \`
72059	kalibera	747	`}`
72142	kalibera	748	`TCROSSPROD_BODY;`
72059	kalibera	749	`}`
		750
72142	kalibera	751	`static void simple_tcrossprod(double *x, int nrx, int ncx,`
		752	`double y, int nry, int ncy, double z)`
		753	`{`
		754	`double sum;`
		755	`TCROSSPROD_BODY;`
		756	`}`
		757
		758
13942	bates	759	`static void matprod(double *x, int nrx, int ncx,`
		760	`double y, int nry, int ncy, double z)`
		761	`{`
60241	ripley	762	`R_xlen_t NRX = nrx, NRY = nry;`
72142	kalibera	763	`if (nrx == 0 \|\| ncx == 0 \|\| nry == 0 \|\| ncy == 0) {`
		764	`/* zero-extent operations should return zeroes */`
		765	`for(R_xlen_t i = 0; i < NRX*ncy; i++) z[i] = 0;`
		766	`return;`
		767	`}`
27297	ripley	768
72142	kalibera	769	`switch(R_Matprod) {`
		770	`case MATPROD_DEFAULT:`
27297	ripley	771	`/* Don't trust the BLAS to handle NA/NaNs correctly: PR#4582`
59068	ripley	772	`* The test is only O(n) here.`
72020	kalibera	773	`*`
		774	`* MKL disclaimer: "LAPACK routines assume that input matrices`
		775	`* do not contain IEEE 754 special values such as INF or NaN values.`
		776	`* Using these special values may cause LAPACK to return unexpected`
		777	`* results or become unstable."`
27297	ripley	778	`*/`
72142	kalibera	779	`if (mayHaveNaNOrInf(x, NRXncx) \|\| mayHaveNaNOrInf(y, NRYncy)) {`
		780	`simple_matprod(x, nrx, ncx, y, nry, ncy, z);`
		781	`return;`
		782	`}`
		783	`break; /* use blas */`
		784	`case MATPROD_INTERNAL:`
		785	`internal_matprod(x, nrx, ncx, y, nry, ncy, z);`
		786	`return;`
		787	`case MATPROD_BLAS:`
		788	`break;`
		789	`case MATPROD_DEFAULT_SIMD:`
		790	`if (mayHaveNaNOrInf_simd(x, NRX*ncx) \|\|`
		791	`mayHaveNaNOrInf_simd(y, NRY*ncy)) {`
		792	`simple_matprod(x, nrx, ncx, y, nry, ncy, z);`
		793	`return;`
		794	`}`
		795	`break; /* use blas */`
		796	`}`
		797
		798	`char transN = "N", transT = "T";`
		799	`double one = 1.0, zero = 0.0;`
		800	`int ione = 1;`
		801
		802	`if (ncy == 1) /* matrix-vector or dot product */`
		803	`F77_CALL(dgemv)(transN, &nrx, &ncx, &one, x,`
76585	ripley	804	`&nrx, y, &ione, &zero, z, &ione FCONE);`
72142	kalibera	805	`else if (nrx == 1) /* vector-matrix */`
		806	`/* Instead of xY, compute (xY)^T == (Y^T)(x^T)`
		807	`The result is a vector, so transposing its content is no-op */`
		808	`F77_CALL(dgemv)(transT, &nry, &ncy, &one, y,`
76585	ripley	809	`&nry, x, &ione, &zero, z, &ione FCONE);`
72142	kalibera	810	`else /* matrix-matrix or outer product */`
76585	ripley	811	`F77_CALL(dgemm)(transN, transN, &nrx, &ncy, &ncx, &one, x,`
		812	`&nrx, y, &nry, &zero, z, &nrx FCONE FCONE);`
27297	ripley	813	`}`
571	ihaka	814
72142	kalibera	815	`static void internal_cmatprod(Rcomplex *x, int nrx, int ncx,`
		816	`Rcomplex y, int nry, int ncy, Rcomplex z)`
2	r	817	`{`
60268	ripley	818	`LDOUBLE sum_i, sum_r;`
72142	kalibera	819	`#define CMATPROD_BODY \`
		820	`int i, j, k; \`
72294	kalibera	821	`double complex xij, yjk; \`
72142	kalibera	822	`R_xlen_t NRX = nrx, NRY = nry; \`
		823	`for (i = 0; i < nrx; i++) \`
		824	`for (k = 0; k < ncy; k++) { \`
		825	`sum_r = 0.0; \`
		826	`sum_i = 0.0; \`
		827	`for (j = 0; j < ncx; j++) { \`
72294	kalibera	828	`xij = toC99(x + (i + j * NRX)); \`
		829	`yjk = toC99(y + (j + k * NRY)); \`
		830	`sum_r += creal(xij * yjk); \`
		831	`sum_i += cimag(xij * yjk); \`
72142	kalibera	832	`} \`
72153	ripley	833	`z[i + k * NRX].r = (double) sum_r; \`
		834	`z[i + k * NRX].i = (double) sum_i; \`
23175	ripley	835	`}`
72142	kalibera	836	`CMATPROD_BODY;`
72083	kalibera	837	`}`
		838
72142	kalibera	839	`static void simple_cmatprod(Rcomplex *x, int nrx, int ncx,`
		840	`Rcomplex y, int nry, int ncy, Rcomplex z)`
72083	kalibera	841	`{`
72142	kalibera	842	`double sum_i, sum_r;`
		843	`CMATPROD_BODY;`
		844	`}`
		845
		846	`static void internal_ccrossprod(Rcomplex *x, int nrx, int ncx,`
		847	`Rcomplex y, int nry, int ncy, Rcomplex z)`
		848	`{`
72083	kalibera	849	`LDOUBLE sum_i, sum_r;`
72142	kalibera	850	`#define CCROSSPROD_BODY \`
		851	`int i, j, k; \`
72294	kalibera	852	`double complex xji, yjk; \`
72142	kalibera	853	`R_xlen_t NRX = nrx, NRY = nry, NCX = ncx; \`
		854	`for (i = 0; i < ncx; i++) \`
		855	`for (k = 0; k < ncy; k++) { \`
		856	`sum_r = 0.0; \`
		857	`sum_i = 0.0; \`
		858	`for (j = 0; j < nrx; j++) { \`
72294	kalibera	859	`xji = toC99(x + (j + i * NRX)); \`
		860	`yjk = toC99(y + (j + k * NRY)); \`
		861	`sum_r += creal(xji * yjk); \`
		862	`sum_i += cimag(xji * yjk); \`
72142	kalibera	863	`} \`
72153	ripley	864	`z[i + k * NCX].r = (double) sum_r; \`
		865	`z[i + k * NCX].i = (double) sum_i; \`
72083	kalibera	866	`}`
72142	kalibera	867	`CCROSSPROD_BODY;`
72083	kalibera	868	`}`
		869
72142	kalibera	870	`static void simple_ccrossprod(Rcomplex *x, int nrx, int ncx,`
		871	`Rcomplex y, int nry, int ncy, Rcomplex z)`
72083	kalibera	872	`{`
72142	kalibera	873	`double sum_i, sum_r;`
		874	`CCROSSPROD_BODY;`
		875	`}`
		876
		877	`static void internal_tccrossprod(Rcomplex *x, int nrx, int ncx,`
		878	`Rcomplex y, int nry, int ncy, Rcomplex z)`
		879	`{`
72083	kalibera	880	`LDOUBLE sum_i, sum_r;`
72142	kalibera	881	`#define TCCROSSPROD_BODY \`
		882	`int i, j, k; \`
72294	kalibera	883	`double complex xij, ykj; \`
72142	kalibera	884	`R_xlen_t NRX = nrx, NRY = nry; \`
		885	`for (i = 0; i < nrx; i++) \`
		886	`for (k = 0; k < nry; k++) { \`
		887	`sum_r = 0.0; \`
		888	`sum_i = 0.0; \`
		889	`for (j = 0; j < ncx; j++) { \`
72294	kalibera	890	`xij = toC99(x + (i + j * NRX)); \`
		891	`ykj = toC99(y + (k + j * NRY)); \`
		892	`sum_r += creal(xij * ykj); \`
		893	`sum_i += cimag(xij * ykj); \`
72142	kalibera	894	`} \`
72153	ripley	895	`z[i + k * NRX].r = (double) sum_r; \`
		896	`z[i + k * NRX].i = (double) sum_i; \`
72083	kalibera	897	`}`
72142	kalibera	898	`TCCROSSPROD_BODY;`
72083	kalibera	899	`}`
		900
72142	kalibera	901	`static void simple_tccrossprod(Rcomplex *x, int nrx, int ncx,`
		902	`Rcomplex y, int nry, int ncy, Rcomplex z)`
		903	`{`
		904	`double sum_i, sum_r;`
		905	`TCCROSSPROD_BODY;`
		906	`}`
		907
72083	kalibera	908	`static void cmatprod(Rcomplex *x, int nrx, int ncx,`
		909	`Rcomplex y, int nry, int ncy, Rcomplex z)`
		910	`{`
		911	`R_xlen_t NRX = nrx, NRY = nry;`
72142	kalibera	912	`if (nrx == 0 \|\| ncx == 0 \|\| nry == 0 \|\| ncy == 0) {`
		913	`/* zero-extent operations should return zeroes */`
		914	`for(R_xlen_t i = 0; i < NRX*ncy; i++) z[i].r = z[i].i = 0;`
		915	`return;`
		916	`}`
72083	kalibera	917
		918	`#ifndef HAVE_FORTRAN_DOUBLE_COMPLEX`
72142	kalibera	919	`if (R_Matprod == MATPROD_INTERNAL)`
		920	`internal_cmatprod(x, nrx, ncx, y, nry, ncy, z);`
		921	`else`
72083	kalibera	922	`simple_cmatprod(x, nrx, ncx, y, nry, ncy, z);`
		923	`#else`
72142	kalibera	924	`switch(R_Matprod) {`
		925	`case MATPROD_DEFAULT:`
		926	`if (cmayHaveNaNOrInf(x, NRXncx) \|\| cmayHaveNaNOrInf(y, NRYncy)) {`
		927	`simple_cmatprod(x, nrx, ncx, y, nry, ncy, z);`
		928	`return;`
		929	`}`
		930	`break; /* use blas */`
		931	`case MATPROD_INTERNAL:`
		932	`internal_cmatprod(x, nrx, ncx, y, nry, ncy, z);`
		933	`return;`
		934	`case MATPROD_BLAS:`
		935	`break;`
		936	`case MATPROD_DEFAULT_SIMD:`
		937	`if (cmayHaveNaNOrInf_simd(x, NRX*ncx) \|\|`
		938	`cmayHaveNaNOrInf_simd(y, NRY*ncy)) {`
		939	`simple_cmatprod(x, nrx, ncx, y, nry, ncy, z);`
		940	`return;`
		941	`}`
		942	`break; /* use blas */`
		943	`}`
		944
		945	`char transa = "N", transb = "N";`
		946	`Rcomplex one, zero;`
		947	`one.r = 1.0; one.i = zero.r = zero.i = 0.0;`
		948
		949	`F77_CALL(zgemm)(transa, transb, &nrx, &ncy, &ncx, &one,`
76585	ripley	950	`x, &nrx, y, &nry, &zero, z, &nrx FCONE FCONE);`
1016	maechler	951	`#endif`
2	r	952	`}`
		953
18753	ripley	954	`static void symcrossprod(double x, int nr, int nc, double z)`
		955	`{`
72142	kalibera	956	`R_xlen_t NR = nr, NC = nc;`
		957	`if (nr == 0 \|\| nc == 0) {`
		958	`/* zero-extent operations should return zeroes */`
60241	ripley	959	`for(R_xlen_t i = 0; i < NC*NC; i++) z[i] = 0;`
72142	kalibera	960	`return;`
18753	ripley	961	`}`
26666	ripley	962
72142	kalibera	963	`switch(R_Matprod) {`
		964	`case MATPROD_DEFAULT:`
		965	`/* see matprod for more details */`
		966	`if (mayHaveNaNOrInf(x, NR*nc)) {`
		967	`simple_crossprod(x, nr, nc, x, nr, nc, z);`
		968	`return;`
		969	`}`
		970	`break; /* use blas */`
		971	`case MATPROD_INTERNAL:`
		972	`internal_crossprod(x, nr, nc, x, nr, nc, z);`
		973	`return;`
		974	`case MATPROD_BLAS:`
		975	`break;`
		976	`case MATPROD_DEFAULT_SIMD:`
		977	`if (mayHaveNaNOrInf_simd(x, NR*nc)) {`
		978	`simple_crossprod(x, nr, nc, x, nr, nc, z);`
		979	`return;`
		980	`}`
		981	`break; /* use blas */`
		982	`}`
		983
		984	`char trans = "T", uplo = "U";`
		985	`double one = 1.0, zero = 0.0;`
		986
76585	ripley	987	`F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc`
		988	`FCONE FCONE);`
72142	kalibera	989	`for (int i = 1; i < nc; i++)`
		990	`for (int j = 0; j < i; j++) z[i + NC j] = z[j + NC i];`
18753	ripley	991	`}`
		992
1820	ihaka	993	`static void crossprod(double *x, int nrx, int ncx,`
		994	`double y, int nry, int ncy, double z)`
2	r	995	`{`
72059	kalibera	996	`R_xlen_t NRX = nrx, NRY = nry;`
72142	kalibera	997	`if (nrx == 0 \|\| ncx == 0 \|\| nry == 0 \|\| ncy == 0) {`
		998	`/* zero-extent operations should return zeroes */`
60241	ripley	999	`R_xlen_t NCX = ncx;`
		1000	`for(R_xlen_t i = 0; i < NCX*ncy; i++) z[i] = 0;`
72142	kalibera	1001	`return;`
26666	ripley	1002	`}`
72142	kalibera	1003
		1004	`switch(R_Matprod) {`
		1005	`case MATPROD_DEFAULT:`
		1006	`/* see matprod for more details */`
		1007	`if (mayHaveNaNOrInf(x, NRXncx) \|\| mayHaveNaNOrInf(y, NRYncy)) {`
		1008	`simple_crossprod(x, nrx, ncx, y, nry, ncy, z);`
		1009	`return;`
		1010	`}`
		1011	`break; /* use blas */`
		1012	`case MATPROD_INTERNAL:`
		1013	`internal_crossprod(x, nrx, ncx, y, nry, ncy, z);`
		1014	`return;`
		1015	`case MATPROD_BLAS:`
		1016	`break;`
		1017	`case MATPROD_DEFAULT_SIMD:`
		1018	`if (mayHaveNaNOrInf_simd(x, NRX*ncx) \|\|`
		1019	`mayHaveNaNOrInf_simd(y, NRY*ncy)) {`
		1020	`simple_crossprod(x, nrx, ncx, y, nry, ncy, z);`
		1021	`return;`
		1022	`}`
		1023	`break; /* use blas */`
		1024	`}`
		1025
		1026	`char transT = "T", transN = "N";`
		1027	`double one = 1.0, zero = 0.0;`
		1028	`int ione = 1;`
		1029
		1030	`if (ncy == 1) /* matrix-vector or dot product */`
		1031	`F77_CALL(dgemv)(transT, &nrx, &ncx, &one, x,`
76585	ripley	1032	`&nrx, y, &ione, &zero, z, &ione FCONE);`
72142	kalibera	1033	`else if (ncx == 1) /* vector-matrix */`
		1034	`/* Instead of (x^T)Y, compute ((x^T)Y)^T == (Y^T)x`
		1035	`The result is a vector, so transposing its content is no-op */`
		1036	`F77_CALL(dgemv)(transT, &nry, &ncy, &one, y,`
76585	ripley	1037	`&nry, x, &ione, &zero, z, &ione FCONE);`
72142	kalibera	1038	`else /* matrix-matrix or outer product */`
		1039	`F77_CALL(dgemm)(transT, transN, &ncx, &ncy, &nrx, &one,`
76585	ripley	1040	`x, &nrx, y, &nry, &zero, z, &ncx FCONE FCONE);`
2	r	1041	`}`
		1042
6994	pd	1043	`static void ccrossprod(Rcomplex *x, int nrx, int ncx,`
		1044	`Rcomplex y, int nry, int ncy, Rcomplex z)`
2	r	1045	`{`
72083	kalibera	1046	`R_xlen_t NRX = nrx, NRY = nry;`
72142	kalibera	1047	`if (nrx == 0 \|\| ncx == 0 \|\| nry == 0 \|\| ncy == 0) {`
		1048	`/* zero-extent operations should return zeroes */`
		1049	`R_xlen_t NCX = ncx;`
		1050	`for(R_xlen_t i = 0; i < NCX*ncy; i++) z[i].r = z[i].i = 0;`
		1051	`return;`
		1052	`}`
23175	ripley	1053
72083	kalibera	1054	`#ifndef HAVE_FORTRAN_DOUBLE_COMPLEX`
72142	kalibera	1055	`if (R_Matprod == MATPROD_INTERNAL)`
		1056	`internal_ccrossprod(x, nrx, ncx, y, nry, ncy, z);`
		1057	`else`
72083	kalibera	1058	`simple_ccrossprod(x, nrx, ncx, y, nry, ncy, z);`
		1059	`#else`
72142	kalibera	1060	`switch(R_Matprod) {`
		1061	`case MATPROD_DEFAULT:`
		1062	`if (cmayHaveNaNOrInf(x, NRXncx) \|\| cmayHaveNaNOrInf(y, NRYncy)) {`
		1063	`simple_ccrossprod(x, nrx, ncx, y, nry, ncy, z);`
		1064	`return;`
		1065	`}`
		1066	`break; /* use blas */`
		1067	`case MATPROD_INTERNAL:`
		1068	`internal_ccrossprod(x, nrx, ncx, y, nry, ncy, z);`
		1069	`return;`
		1070	`case MATPROD_BLAS:`
		1071	`break;`
		1072	`case MATPROD_DEFAULT_SIMD:`
		1073	`if (cmayHaveNaNOrInf_simd(x, NRX*ncx) \|\|`
		1074	`cmayHaveNaNOrInf_simd(y, NRY*ncy)) {`
		1075	`simple_ccrossprod(x, nrx, ncx, y, nry, ncy, z);`
		1076	`return;`
		1077	`}`
		1078	`break; /* use blas */`
		1079	`}`
		1080
		1081	`char transa = "T", transb = "N";`
		1082	`Rcomplex one, zero;`
		1083	`one.r = 1.0; one.i = zero.r = zero.i = 0.0;`
		1084
		1085	`F77_CALL(zgemm)(transa, transb, &ncx, &ncy, &nrx, &one,`
76585	ripley	1086	`x, &nrx, y, &nry, &zero, z, &ncx FCONE FCONE);`
72083	kalibera	1087	`#endif`
2	r	1088	`}`
36679	ripley	1089
		1090	`static void symtcrossprod(double x, int nr, int nc, double z)`
		1091	`{`
72142	kalibera	1092	`R_xlen_t NR = nr;`
		1093	`if (nr == 0 \|\| nc == 0) {`
		1094	`/* zero-extent operations should return zeroes */`
60241	ripley	1095	`for(R_xlen_t i = 0; i < NR*NR; i++) z[i] = 0;`
72142	kalibera	1096	`return;`
36679	ripley	1097	`}`
		1098
72142	kalibera	1099	`switch(R_Matprod) {`
		1100	`case MATPROD_DEFAULT:`
		1101	`/* see matprod for more details */`
		1102	`if (mayHaveNaNOrInf(x, NR*nc)) {`
		1103	`simple_tcrossprod(x, nr, nc, x, nr, nc, z);`
		1104	`return;`
		1105	`}`
		1106	`break; /* use blas */`
		1107	`case MATPROD_INTERNAL:`
		1108	`internal_tcrossprod(x, nr, nc, x, nr, nc, z);`
		1109	`return;`
		1110	`case MATPROD_BLAS:`
		1111	`break;`
		1112	`case MATPROD_DEFAULT_SIMD:`
		1113	`if (mayHaveNaNOrInf_simd(x, NR*nc)) {`
		1114	`simple_tcrossprod(x, nr, nc, x, nr, nc, z);`
		1115	`return;`
		1116	`}`
		1117	`break; /* use blas */`
		1118	`}`
		1119
		1120	`char trans = "N", uplo = "U";`
		1121	`double one = 1.0, zero = 0.0;`
		1122
76585	ripley	1123	`F77_CALL(dsyrk)(uplo, trans, &nr, &nc, &one, x, &nr, &zero, z, &nr`
		1124	`FCONE FCONE);`
72142	kalibera	1125	`for (int i = 1; i < nr; i++)`
		1126	`for (int j = 0; j < i; j++) z[i + nr j] = z[j + nr i];`
36679	ripley	1127	`}`
		1128
		1129	`static void tcrossprod(double *x, int nrx, int ncx,`
		1130	`double y, int nry, int ncy, double z)`
		1131	`{`
72142	kalibera	1132	`R_xlen_t NRX = nrx, NRY = nry;`
		1133	`if (nrx == 0 \|\| ncx == 0 \|\| nry == 0 \|\| ncy == 0) {`
		1134	`/* zero-extent operations should return zeroes */`
		1135	`for(R_xlen_t i = 0; i < NRX*nry; i++) z[i] = 0;`
		1136	`return;`
		1137	`}`
		1138
		1139	`switch(R_Matprod) {`
		1140	`case MATPROD_DEFAULT:`
		1141	`if (mayHaveNaNOrInf(x, NRXncx) \|\| mayHaveNaNOrInf(y, NRYncy)) {`
		1142	`simple_tcrossprod(x, nrx, ncx, y, nry, ncy, z);`
		1143	`return;`
		1144	`}`
		1145	`break; /* use blas */`
		1146	`case MATPROD_INTERNAL:`
		1147	`internal_tcrossprod(x, nrx, ncx, y, nry, ncy, z);`
		1148	`return;`
		1149	`case MATPROD_BLAS:`
		1150	`break;`
		1151	`case MATPROD_DEFAULT_SIMD:`
		1152	`if (mayHaveNaNOrInf_simd(x, NRX*ncx) \|\|`
		1153	`mayHaveNaNOrInf_simd(y, NRY*ncy)) {`
		1154	`simple_tcrossprod(x, nrx, ncx, y, nry, ncy, z);`
		1155	`return;`
		1156	`}`
		1157	`break; /* use blas */`
		1158	`}`
		1159
72036	kalibera	1160	`char transN = "N", transT = "T";`
36679	ripley	1161	`double one = 1.0, zero = 0.0;`
72036	kalibera	1162	`int ione = 1;`
72142	kalibera	1163
		1164	`if (nry == 1) /* matrix-vector or dot product */`
		1165	`F77_CALL(dgemv)(transN, &nrx, &ncx, &one, x,`
76585	ripley	1166	`&nrx, y, &ione, &zero, z, &ione FCONE);`
72142	kalibera	1167	`else if (nrx == 1) /* vector-matrix */`
		1168	`/* Instead of x(Y^T), compute (x(Y^T))^T == Y(x^T)`
		1169	`The result is a vector, so transposing its content is no-op */`
		1170	`F77_CALL(dgemv)(transN, &nry, &ncy, &one, y,`
76585	ripley	1171	`&nry, x, &ione, &zero, z, &ione FCONE);`
72142	kalibera	1172	`else /* matrix-matrix or outer product */`
		1173	`F77_CALL(dgemm)(transN, transT, &nrx, &nry, &ncx, &one,`
76585	ripley	1174	`x, &nrx, y, &nry, &zero, z, &nrx FCONE FCONE);`
36679	ripley	1175	`}`
		1176
		1177	`static void tccrossprod(Rcomplex *x, int nrx, int ncx,`
		1178	`Rcomplex y, int nry, int ncy, Rcomplex z)`
		1179	`{`
72083	kalibera	1180	`R_xlen_t NRX = nrx, NRY = nry;`
72142	kalibera	1181	`if (nrx == 0 \|\| ncx == 0 \|\| nry == 0 \|\| ncy == 0) {`
		1182	`/* zero-extent operations should return zeroes */`
		1183	`for(R_xlen_t i = 0; i < NRX*nry; i++) z[i].r = z[i].i = 0;`
		1184	`return;`
		1185	`}`
36679	ripley	1186
72083	kalibera	1187	`#ifndef HAVE_FORTRAN_DOUBLE_COMPLEX`
72142	kalibera	1188	`if (R_Matprod == MATPROD_INTERNAL)`
		1189	`internal_tccrossprod(x, nrx, ncx, y, nry, ncy, z);`
		1190	`else`
72083	kalibera	1191	`simple_tccrossprod(x, nrx, ncx, y, nry, ncy, z);`
		1192	`#else`
72142	kalibera	1193	`switch(R_Matprod) {`
		1194	`case MATPROD_DEFAULT:`
		1195	`if (cmayHaveNaNOrInf(x, NRXncx) \|\| cmayHaveNaNOrInf(y, NRYncy)) {`
		1196	`simple_tccrossprod(x, nrx, ncx, y, nry, ncy, z);`
		1197	`return;`
		1198	`}`
		1199	`break; /* use blas */`
		1200	`case MATPROD_INTERNAL:`
		1201	`internal_tccrossprod(x, nrx, ncx, y, nry, ncy, z);`
		1202	`return;`
		1203	`case MATPROD_BLAS:`
		1204	`break;`
		1205	`case MATPROD_DEFAULT_SIMD:`
		1206	`if (cmayHaveNaNOrInf_simd(x, NRX*ncx) \|\|`
		1207	`cmayHaveNaNOrInf_simd(y, NRY*ncy)) {`
		1208	`simple_tccrossprod(x, nrx, ncx, y, nry, ncy, z);`
		1209	`return;`
		1210	`}`
		1211	`break; /* use blas */`
		1212	`}`
		1213
		1214	`char transa = "N", transb = "T";`
		1215	`Rcomplex one, zero;`
		1216	`one.r = 1.0; one.i = zero.r = zero.i = 0.0;`
		1217
		1218	`F77_CALL(zgemm)(transa, transb, &nrx, &nry, &ncx, &one,`
76585	ripley	1219	`x, &nrx, y, &nry, &zero, z, &nrx FCONE FCONE);`
72083	kalibera	1220	`#endif`
36679	ripley	1221	`}`
		1222
		1223
		1224	`/* "%%" (op = 0), crossprod (op = 1) or tcrossprod (op = 2) /`
36990	ripley	1225	`SEXP attribute_hidden do_matprod(SEXP call, SEXP op, SEXP args, SEXP rho)`
2	r	1226	`{`
1820	ihaka	1227	`int ldx, ldy, nrx, ncx, nry, ncy, mode;`
18753	ripley	1228	`SEXP x = CAR(args), y = CADR(args), xdims, ydims, ans;`
		1229	`Rboolean sym;`
2	r	1230
59239	ripley	1231	`if (PRIMVAL(op) == 0 && /* %% is primitive, the others are .Internal() /`
45446	ripley	1232	`(IS_S4_OBJECT(x) \|\| IS_S4_OBJECT(y))`
41281	ripley	1233	`&& R_has_methods(op)) {`
41365	ripley	1234	`SEXP s, value;`
		1235	`/* Remove argument names to ensure positional matching */`
		1236	`for(s = args; s != R_NilValue; s = CDR(s)) SET_TAG(s, R_NilValue);`
42292	ripley	1237	`value = R_possible_dispatch(call, op, args, rho, FALSE);`
59239	ripley	1238	`if (value) return value;`
21875	jmc	1239	`}`
		1240
69326	luke	1241	`checkArity(op, args);`
18753	ripley	1242	`sym = isNull(y);`
36679	ripley	1243	`if (sym && (PRIMVAL(op) > 0)) y = x;`
18753	ripley	1244	`if ( !(isNumeric(x) \|\| isComplex(x)) \|\| !(isNumeric(y) \|\| isComplex(y)) )`
48092	maechler	1245	`errorcall(call, _("requires numeric/complex matrix/vector arguments"));`
2	r	1246
1820	ihaka	1247	`xdims = getAttrib(x, R_DimSymbol);`
		1248	`ydims = getAttrib(y, R_DimSymbol);`
		1249	`ldx = length(xdims);`
		1250	`ldy = length(ydims);`
2	r	1251
26472	ihaka	1252	`if (ldx != 2 && ldy != 2) { /* x and y non-matrices */`
66495	maechler	1253	`// for crossprod, allow two cases: n x n ==> (1,n) x (n,1); 1 x n = (n, 1) x (1, n)`
		1254	`if (PRIMVAL(op) == 1 && LENGTH(x) == 1) {`
		1255	`nrx = ncx = nry = 1;`
		1256	`ncy = LENGTH(y);`
2	r	1257	`}`
1820	ihaka	1258	`else {`
66495	maechler	1259	`nry = LENGTH(y);`
		1260	`ncy = 1;`
		1261	`if (PRIMVAL(op) == 0) {`
		1262	`nrx = 1;`
		1263	`ncx = LENGTH(x);`
		1264	`if(ncx == 1) { // y as row vector`
		1265	`ncy = nry;`
		1266	`nry = 1;`
		1267	`}`
		1268	`}`
		1269	`else {`
		1270	`nrx = LENGTH(x);`
		1271	`ncx = 1;`
		1272	`}`
2	r	1273	`}`
1820	ihaka	1274	`}`
26472	ihaka	1275	`else if (ldx != 2) { /* x not a matrix */`
1820	ihaka	1276	`nry = INTEGER(ydims)[0];`
		1277	`ncy = INTEGER(ydims)[1];`
		1278	`nrx = 0;`
		1279	`ncx = 0;`
1839	ihaka	1280	`if (PRIMVAL(op) == 0) {`
26472	ihaka	1281	`if (LENGTH(x) == nry) { /* x as row vector */`
1820	ihaka	1282	`nrx = 1;`
48092	maechler	1283	`ncx = nry; /* == LENGTH(x) */`
2213	maechler	1284	`}`
26472	ihaka	1285	`else if (nry == 1) { /* x as col vector */`
1820	ihaka	1286	`nrx = LENGTH(x);`
		1287	`ncx = 1;`
		1288	`}`
2	r	1289	`}`
48574	maechler	1290	`else if (PRIMVAL(op) == 1) { /* crossprod() */`
36679	ripley	1291	`if (LENGTH(x) == nry) { /* x is a col vector */`
48092	maechler	1292	`nrx = nry; /* == LENGTH(x) */`
1820	ihaka	1293	`ncx = 1;`
		1294	`}`
50989	maechler	1295	`/* else if (nry == 1) ... not being too tolerant`
		1296	`to treat x as row vector, as t(x) is row vector */`
2	r	1297	`}`
48574	maechler	1298	`else { /* tcrossprod */`
		1299	`if (LENGTH(x) == ncy) { /* x as row vector */`
		1300	`nrx = 1;`
		1301	`ncx = ncy; /* == LENGTH(x) */`
		1302	`}`
50989	maechler	1303	`else if (ncy == 1) { /* x as col vector */`
		1304	`nrx = LENGTH(x);`
		1305	`ncx = 1;`
		1306	`}`
48574	maechler	1307	`}`
1820	ihaka	1308	`}`
26472	ihaka	1309	`else if (ldy != 2) { /* y not a matrix */`
1820	ihaka	1310	`nrx = INTEGER(xdims)[0];`
		1311	`ncx = INTEGER(xdims)[1];`
		1312	`nry = 0;`
		1313	`ncy = 0;`
1839	ihaka	1314	`if (PRIMVAL(op) == 0) {`
26472	ihaka	1315	`if (LENGTH(y) == ncx) { /* y as col vector */`
50989	maechler	1316	`nry = ncx;`
1820	ihaka	1317	`ncy = 1;`
		1318	`}`
26472	ihaka	1319	`else if (ncx == 1) { /* y as row vector */`
1820	ihaka	1320	`nry = 1;`
		1321	`ncy = LENGTH(y);`
		1322	`}`
2	r	1323	`}`
50989	maechler	1324	`else if (PRIMVAL(op) == 1) { /* crossprod() */`
36679	ripley	1325	`if (LENGTH(y) == nrx) { /* y is a col vector */`
50989	maechler	1326	`nry = nrx;`
1820	ihaka	1327	`ncy = 1;`
66495	maechler	1328	`} else if (nrx == 1) { // y as row vector`
		1329	`nry = 1;`
		1330	`ncy = LENGTH(y);`
1820	ihaka	1331	`}`
2	r	1332	`}`
66495	maechler	1333	`else { // tcrossprod`
		1334	`if (nrx == 1) { // y as row vector`
		1335	`nry = 1;`
		1336	`ncy = LENGTH(y);`
		1337	`}`
		1338	`else { // y is a col vector`
		1339	`nry = LENGTH(y);`
		1340	`ncy = 1;`
		1341	`}`
50989	maechler	1342	`}`
1820	ihaka	1343	`}`
26472	ihaka	1344	`else { /* x and y matrices */`
1820	ihaka	1345	`nrx = INTEGER(xdims)[0];`
		1346	`ncx = INTEGER(xdims)[1];`
		1347	`nry = INTEGER(ydims)[0];`
		1348	`ncy = INTEGER(ydims)[1];`
		1349	`}`
32350	maechler	1350	`/* nr[ow](.) and nc[ol](.) are now defined for x and y */`
2	r	1351
1839	ihaka	1352	`if (PRIMVAL(op) == 0) {`
41686	ripley	1353	`/* primitive, so use call */`
1839	ihaka	1354	`if (ncx != nry)`
32860	ripley	1355	`errorcall(call, _("non-conformable arguments"));`
1820	ihaka	1356	`}`
36679	ripley	1357	`else if (PRIMVAL(op) == 1) {`
1839	ihaka	1358	`if (nrx != nry)`
41686	ripley	1359	`error(_("non-conformable arguments"));`
1820	ihaka	1360	`}`
36679	ripley	1361	`else {`
		1362	`if (ncx != ncy)`
41686	ripley	1363	`error(_("non-conformable arguments"));`
36679	ripley	1364	`}`
1820	ihaka	1365
1839	ihaka	1366	`if (isComplex(CAR(args)) \|\| isComplex(CADR(args)))`
1820	ihaka	1367	`mode = CPLXSXP;`
		1368	`else`
		1369	`mode = REALSXP;`
10172	luke	1370	`SETCAR(args, coerceVector(CAR(args), mode));`
		1371	`SETCADR(args, coerceVector(CADR(args), mode));`
1820	ihaka	1372
32350	maechler	1373	`if (PRIMVAL(op) == 0) { /* op == 0 : matprod() */`
26472	ihaka	1374
1820	ihaka	1375	`PROTECT(ans = allocMatrix(mode, nrx, ncy));`
1839	ihaka	1376	`if (mode == CPLXSXP)`
1820	ihaka	1377	`cmatprod(COMPLEX(CAR(args)), nrx, ncx,`
		1378	`COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));`
2	r	1379	`else`
1820	ihaka	1380	`matprod(REAL(CAR(args)), nrx, ncx,`
		1381	`REAL(CADR(args)), nry, ncy, REAL(ans));`
26472	ihaka	1382
1820	ihaka	1383	`PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));`
		1384	`PROTECT(ydims = getAttrib(CADR(args), R_DimNamesSymbol));`
26472	ihaka	1385
1820	ihaka	1386	`if (xdims != R_NilValue \|\| ydims != R_NilValue) {`
32350	maechler	1387	`SEXP dimnames, dimnamesnames, dnx=R_NilValue, dny=R_NilValue;`
		1388
		1389	`/* allocate dimnames and dimnamesnames */`
		1390
6994	pd	1391	`PROTECT(dimnames = allocVector(VECSXP, 2));`
		1392	`PROTECT(dimnamesnames = allocVector(STRSXP, 2));`
		1393	`if (xdims != R_NilValue) {`
32350	maechler	1394	`if (ldx == 2 \|\| ncx == 1) {`
10172	luke	1395	`SET_VECTOR_ELT(dimnames, 0, VECTOR_ELT(xdims, 0));`
32350	maechler	1396	`dnx = getAttrib(xdims, R_NamesSymbol);`
		1397	`if(!isNull(dnx))`
		1398	`SET_STRING_ELT(dimnamesnames, 0, STRING_ELT(dnx, 0));`
8628	pd	1399	`}`
6994	pd	1400	`}`
32350	maechler	1401
32356	maechler	1402	`#define YDIMS_ET_CETERA \`
		1403	`if (ydims != R_NilValue) { \`
		1404	`if (ldy == 2) { \`
		1405	`SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 1)); \`
		1406	`dny = getAttrib(ydims, R_NamesSymbol); \`
		1407	`if(!isNull(dny)) \`
		1408	`SET_STRING_ELT(dimnamesnames, 1, STRING_ELT(dny, 1)); \`
		1409	`} else if (nry == 1) { \`
		1410	`SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 0)); \`
		1411	`dny = getAttrib(ydims, R_NamesSymbol); \`
		1412	`if(!isNull(dny)) \`
		1413	`SET_STRING_ELT(dimnamesnames, 1, STRING_ELT(dny, 0)); \`
		1414	`} \`
		1415	`} \`
		1416	`\`
		1417	`/* We sometimes attach a dimnames attribute \`
		1418	`* whose elements are all NULL ... \`
		1419	`* This is ugly but causes no real damage. \`
		1420	`* Now (2.1.0 ff), we don't anymore: */ \`
		1421	`if (VECTOR_ELT(dimnames,0) != R_NilValue \|\| \`
		1422	`VECTOR_ELT(dimnames,1) != R_NilValue) { \`
		1423	`if (dnx != R_NilValue \|\| dny != R_NilValue) \`
		1424	`setAttrib(dimnames, R_NamesSymbol, dimnamesnames); \`
		1425	`setAttrib(ans, R_DimNamesSymbol, dimnames); \`
		1426	`} \`
32350	maechler	1427	`UNPROTECT(2)`
		1428
		1429	`YDIMS_ET_CETERA;`
2	r	1430	`}`
1820	ihaka	1431	`}`
26472	ihaka	1432
36679	ripley	1433	`else if (PRIMVAL(op) == 1) { /* op == 1: crossprod() */`
26472	ihaka	1434
1820	ihaka	1435	`PROTECT(ans = allocMatrix(mode, ncx, ncy));`
1839	ihaka	1436	`if (mode == CPLXSXP)`
23175	ripley	1437	`if(sym)`
		1438	`ccrossprod(COMPLEX(CAR(args)), nrx, ncx,`
		1439	`COMPLEX(CAR(args)), nry, ncy, COMPLEX(ans));`
		1440	`else`
		1441	`ccrossprod(COMPLEX(CAR(args)), nrx, ncx,`
		1442	`COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));`
18753	ripley	1443	`else {`
		1444	`if(sym)`
		1445	`symcrossprod(REAL(CAR(args)), nrx, ncx, REAL(ans));`
		1446	`else`
29387	ripley	1447	`crossprod(REAL(CAR(args)), nrx, ncx,`
		1448	`REAL(CADR(args)), nry, ncy, REAL(ans));`
18753	ripley	1449	`}`
26472	ihaka	1450
1820	ihaka	1451	`PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));`
26472	ihaka	1452	`if (sym)`
18775	ripley	1453	`PROTECT(ydims = xdims);`
		1454	`else`
		1455	`PROTECT(ydims = getAttrib(CADR(args), R_DimNamesSymbol));`
26472	ihaka	1456
1820	ihaka	1457	`if (xdims != R_NilValue \|\| ydims != R_NilValue) {`
7081	pd	1458	`SEXP dimnames, dimnamesnames, dnx=R_NilValue, dny=R_NilValue;`
26472	ihaka	1459
		1460	`/* allocate dimnames and dimnamesnames */`
		1461
6994	pd	1462	`PROTECT(dimnames = allocVector(VECSXP, 2));`
		1463	`PROTECT(dimnamesnames = allocVector(STRSXP, 2));`
26472	ihaka	1464
6994	pd	1465	`if (xdims != R_NilValue) {`
32350	maechler	1466	`if (ldx == 2) {/* not nrx==1 : .. fixed, ihaka 2003-09-30 */`
		1467	`SET_VECTOR_ELT(dimnames, 0, VECTOR_ELT(xdims, 1));`
26472	ihaka	1468	`dnx = getAttrib(xdims, R_NamesSymbol);`
		1469	`if(!isNull(dnx))`
		1470	`SET_STRING_ELT(dimnamesnames, 0, STRING_ELT(dnx, 1));`
		1471	`}`
6994	pd	1472	`}`
26472	ihaka	1473
32350	maechler	1474	`YDIMS_ET_CETERA;`
2	r	1475	`}`
36679	ripley	1476
1820	ihaka	1477	`}`
36679	ripley	1478	`else { /* op == 2: tcrossprod() */`
		1479
		1480	`PROTECT(ans = allocMatrix(mode, nrx, nry));`
		1481	`if (mode == CPLXSXP)`
		1482	`if(sym)`
		1483	`tccrossprod(COMPLEX(CAR(args)), nrx, ncx,`
		1484	`COMPLEX(CAR(args)), nry, ncy, COMPLEX(ans));`
		1485	`else`
		1486	`tccrossprod(COMPLEX(CAR(args)), nrx, ncx,`
		1487	`COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));`
		1488	`else {`
		1489	`if(sym)`
		1490	`symtcrossprod(REAL(CAR(args)), nrx, ncx, REAL(ans));`
		1491	`else`
		1492	`tcrossprod(REAL(CAR(args)), nrx, ncx,`
		1493	`REAL(CADR(args)), nry, ncy, REAL(ans));`
		1494	`}`
		1495
		1496	`PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));`
		1497	`if (sym)`
		1498	`PROTECT(ydims = xdims);`
		1499	`else`
		1500	`PROTECT(ydims = getAttrib(CADR(args), R_DimNamesSymbol));`
		1501
		1502	`if (xdims != R_NilValue \|\| ydims != R_NilValue) {`
		1503	`SEXP dimnames, dimnamesnames, dnx=R_NilValue, dny=R_NilValue;`
		1504
		1505	`/* allocate dimnames and dimnamesnames */`
		1506
		1507	`PROTECT(dimnames = allocVector(VECSXP, 2));`
		1508	`PROTECT(dimnamesnames = allocVector(STRSXP, 2));`
		1509
		1510	`if (xdims != R_NilValue) {`
		1511	`if (ldx == 2) {`
		1512	`SET_VECTOR_ELT(dimnames, 0, VECTOR_ELT(xdims, 0));`
		1513	`dnx = getAttrib(xdims, R_NamesSymbol);`
		1514	`if(!isNull(dnx))`
		1515	`SET_STRING_ELT(dimnamesnames, 0, STRING_ELT(dnx, 0));`
		1516	`}`
		1517	`}`
		1518	`if (ydims != R_NilValue) {`
		1519	`if (ldy == 2) {`
		1520	`SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 0));`
		1521	`dny = getAttrib(ydims, R_NamesSymbol);`
		1522	`if(!isNull(dny))`
		1523	`SET_STRING_ELT(dimnamesnames, 1, STRING_ELT(dny, 0));`
		1524	`}`
		1525	`}`
		1526	`if (VECTOR_ELT(dimnames,0) != R_NilValue \|\|`
		1527	`VECTOR_ELT(dimnames,1) != R_NilValue) {`
		1528	`if (dnx != R_NilValue \|\| dny != R_NilValue)`
		1529	`setAttrib(dimnames, R_NamesSymbol, dimnamesnames);`
		1530	`setAttrib(ans, R_DimNamesSymbol, dimnames);`
		1531	`}`
		1532
		1533	`UNPROTECT(2);`
		1534	`}`
		1535	`}`
1820	ihaka	1536	`UNPROTECT(3);`
		1537	`return ans;`
2	r	1538	`}`
32350	maechler	1539	`#undef YDIMS_ET_CETERA`
2	r	1540
36990	ripley	1541	`SEXP attribute_hidden do_transpose(SEXP call, SEXP op, SEXP args, SEXP rho)`
2	r	1542	`{`
59052	ripley	1543	`SEXP a, r, dims, dimnames, dimnamesnames = R_NilValue,`
6994	pd	1544	`ndimnamesnames, rnames, cnames;`
59052	ripley	1545	`int ldim, ncol = 0, nrow = 0;`
		1546	`R_xlen_t len = 0;`
2	r	1547
1820	ihaka	1548	`checkArity(op, args);`
		1549	`a = CAR(args);`
2	r	1550
1839	ihaka	1551	`if (isVector(a)) {`
1820	ihaka	1552	`dims = getAttrib(a, R_DimSymbol);`
32350	maechler	1553	`ldim = length(dims);`
1858	ihaka	1554	`rnames = R_NilValue;`
		1555	`cnames = R_NilValue;`
32350	maechler	1556	`switch(ldim) {`
1820	ihaka	1557	`case 0:`
59086	ripley	1558	`len = nrow = LENGTH(a);`
1858	ihaka	1559	`ncol = 1;`
		1560	`rnames = getAttrib(a, R_NamesSymbol);`
36597	maechler	1561	`dimnames = rnames;/* for isNull() below*/`
1858	ihaka	1562	`break;`
1820	ihaka	1563	`case 1:`
59086	ripley	1564	`len = nrow = LENGTH(a);`
1820	ihaka	1565	`ncol = 1;`
32350	maechler	1566	`dimnames = getAttrib(a, R_DimNamesSymbol);`
		1567	`if (dimnames != R_NilValue) {`
		1568	`rnames = VECTOR_ELT(dimnames, 0);`
		1569	`dimnamesnames = getAttrib(dimnames, R_NamesSymbol);`
		1570	`}`
1820	ihaka	1571	`break;`
		1572	`case 2:`
		1573	`ncol = ncols(a);`
		1574	`nrow = nrows(a);`
59052	ripley	1575	`len = XLENGTH(a);`
1858	ihaka	1576	`dimnames = getAttrib(a, R_DimNamesSymbol);`
		1577	`if (dimnames != R_NilValue) {`
10172	luke	1578	`rnames = VECTOR_ELT(dimnames, 0);`
		1579	`cnames = VECTOR_ELT(dimnames, 1);`
6994	pd	1580	`dimnamesnames = getAttrib(dimnames, R_NamesSymbol);`
1858	ihaka	1581	`}`
1820	ihaka	1582	`break;`
		1583	`default:`
		1584	`goto not_matrix;`
2	r	1585	`}`
1820	ihaka	1586	`}`
1858	ihaka	1587	`else`
		1588	`goto not_matrix;`
68092	luke	1589	`PROTECT(dimnamesnames);`
1820	ihaka	1590	`PROTECT(r = allocVector(TYPEOF(a), len));`
59052	ripley	1591	`R_xlen_t i, j, l_1 = len-1;`
1820	ihaka	1592	`switch (TYPEOF(a)) {`
		1593	`case LGLSXP:`
		1594	`case INTSXP:`
52876	maechler	1595	`// filling in columnwise, "accessing row-wise":`
68923	ripley	1596	`for (i = 0, j = 0; i < len; i++, j += nrow) {`
		1597	`if (j > l_1) j -= l_1;`
		1598	`INTEGER(r)[i] = INTEGER(a)[j];`
		1599	`}`
		1600	`break;`
1820	ihaka	1601	`case REALSXP:`
68923	ripley	1602	`for (i = 0, j = 0; i < len; i++, j += nrow) {`
		1603	`if (j > l_1) j -= l_1;`
		1604	`REAL(r)[i] = REAL(a)[j];`
		1605	`}`
		1606	`break;`
1820	ihaka	1607	`case CPLXSXP:`
68923	ripley	1608	`for (i = 0, j = 0; i < len; i++, j += nrow) {`
		1609	`if (j > l_1) j -= l_1;`
		1610	`COMPLEX(r)[i] = COMPLEX(a)[j];`
		1611	`}`
		1612	`break;`
1820	ihaka	1613	`case STRSXP:`
68923	ripley	1614	`for (i = 0, j = 0; i < len; i++, j += nrow) {`
		1615	`if (j > l_1) j -= l_1;`
		1616	`SET_STRING_ELT(r, i, STRING_ELT(a,j));`
		1617	`}`
		1618	`break;`
1858	ihaka	1619	`case VECSXP:`
68923	ripley	1620	`for (i = 0, j = 0; i < len; i++, j += nrow) {`
		1621	`if (j > l_1) j -= l_1;`
		1622	`SET_VECTOR_ELT(r, i, VECTOR_ELT(a,j));`
		1623	`}`
		1624	`break;`
29902	ripley	1625	`case RAWSXP:`
68923	ripley	1626	`for (i = 0, j = 0; i < len; i++, j += nrow) {`
		1627	`if (j > l_1) j -= l_1;`
		1628	`RAW(r)[i] = RAW(a)[j];`
		1629	`}`
		1630	`break;`
1858	ihaka	1631	`default:`
68923	ripley	1632	`UNPROTECT(2); /* r, dimnamesnames */`
		1633	`goto not_matrix;`
1820	ihaka	1634	`}`
1858	ihaka	1635	`PROTECT(dims = allocVector(INTSXP, 2));`
1820	ihaka	1636	`INTEGER(dims)[0] = ncol;`
		1637	`INTEGER(dims)[1] = nrow;`
		1638	`setAttrib(r, R_DimSymbol, dims);`
68092	luke	1639	`UNPROTECT(1); /* dims */`
36597	maechler	1640	`/* R <= 2.2.0: dropped list(NULL,NULL) dimnames :`
		1641	`* if(rnames != R_NilValue \|\| cnames != R_NilValue) */`
		1642	`if(!isNull(dimnames)) {`
1858	ihaka	1643	`PROTECT(dimnames = allocVector(VECSXP, 2));`
10172	luke	1644	`SET_VECTOR_ELT(dimnames, 0, cnames);`
		1645	`SET_VECTOR_ELT(dimnames, 1, rnames);`
6994	pd	1646	`if(!isNull(dimnamesnames)) {`
		1647	`PROTECT(ndimnamesnames = allocVector(VECSXP, 2));`
38530	ripley	1648	`SET_VECTOR_ELT(ndimnamesnames, 1, STRING_ELT(dimnamesnames, 0));`
		1649	`SET_VECTOR_ELT(ndimnamesnames, 0,`
32350	maechler	1650	`(ldim == 2) ? STRING_ELT(dimnamesnames, 1):`
		1651	`R_BlankString);`
6994	pd	1652	`setAttrib(dimnames, R_NamesSymbol, ndimnamesnames);`
68092	luke	1653	`UNPROTECT(1); /* ndimnamesnames */`
6994	pd	1654	`}`
1858	ihaka	1655	`setAttrib(r, R_DimNamesSymbol, dimnames);`
68092	luke	1656	`UNPROTECT(1); /* dimnames */`
1820	ihaka	1657	`}`
1858	ihaka	1658	`copyMostAttrib(a, r);`
68092	luke	1659	`UNPROTECT(2); /* r, dimnamesnames */`
1820	ihaka	1660	`return r;`
		1661	`not_matrix:`
41686	ripley	1662	`error(_("argument is not a matrix"));`
1820	ihaka	1663	`return call;/* never used; just for -Wall */`
2	r	1664	`}`
		1665
13214	maechler	1666	`/*`
		1667	`New version of aperm, using strides for speed.`
		1668	`Jonathan Rougier <J.C.Rougier@durham.ac.uk>`
		1669
		1670	`v1.0 30.01.01`
		1671
		1672	`M.Maechler : expanded all ../include/Rdefines.h macros`
		1673	`*/`
		1674
		1675	`/* this increments iip and sets j using strides */`
		1676
		1677	`#define CLICKJ \`
52649	ripley	1678	`for (itmp = 0; itmp < n; itmp++) \`
		1679	`if (iip[itmp] == isr[itmp]-1) iip[itmp] = 0; \`
13214	maechler	1680	`else { \`
		1681	`iip[itmp]++; \`
		1682	`break; \`
		1683	`} \`
68923	ripley	1684	`for (lj = 0, itmp = 0; itmp < n; itmp++) \`
60305	ripley	1685	`lj += iip[itmp] * stride[itmp];`
17515	ripley	1686
13214	maechler	1687	`/* aperm (a, perm, resize = TRUE) */`
36990	ripley	1688	`SEXP attribute_hidden do_aperm(SEXP call, SEXP op, SEXP args, SEXP rho)`
2	r	1689	`{`
52649	ripley	1690	`SEXP a, perm, r, dimsa, dimsr, dna;`
60305	ripley	1691	`int i, j, n, itmp;`
2	r	1692
1820	ihaka	1693	`checkArity(op, args);`
2	r	1694
1820	ihaka	1695	`a = CAR(args);`
13214	maechler	1696	`if (!isArray(a))`
41686	ripley	1697	`error(_("invalid first argument, must be an array"));`
13214	maechler	1698
1820	ihaka	1699	`PROTECT(dimsa = getAttrib(a, R_DimSymbol));`
13214	maechler	1700	`n = LENGTH(dimsa);`
52649	ripley	1701	`int *isa = INTEGER(dimsa);`
2	r	1702
13214	maechler	1703	`/* check the permutation */`
		1704
55420	luke	1705	`int pp = (int ) R_alloc((size_t) n, sizeof(int));`
52649	ripley	1706	`perm = CADR(args);`
13214	maechler	1707	`if (length(perm) == 0) {`
52649	ripley	1708	`for (i = 0; i < n; i++) pp[i] = n-1-i;`
		1709	`} else {`
61696	maechler	1710	`if (LENGTH(perm) != n)`
		1711	`error(_("'perm' is of wrong length %d (!= %d)"),`
		1712	`LENGTH(perm), n);`
		1713	`if (isString(perm)) {`
		1714	`SEXP dna = getAttrib(a, R_DimNamesSymbol);`
		1715	`if (isNull(dna))`
		1716	`error(_("'a' does not have named dimnames"));`
		1717	`SEXP dnna = getAttrib(dna, R_NamesSymbol);`
		1718	`if (isNull(dnna))`
		1719	`error(_("'a' does not have named dimnames"));`
		1720	`for (i = 0; i < n; i++) {`
		1721	`const char *this = translateChar(STRING_ELT(perm, i));`
		1722	`for (j = 0; j < n; j++)`
		1723	`if (streql(translateChar(STRING_ELT(dnna, j)),`
		1724	`this)) {pp[i] = j; break;}`
		1725	`if (j >= n)`
		1726	`error(_("'perm[%d]' does not match a dimension name"), i+1);`
		1727	`}`
		1728	`} else {`
		1729	`PROTECT(perm = coerceVector(perm, INTSXP));`
52649	ripley	1730	`for (i = 0; i < n; i++) pp[i] = INTEGER(perm)[i] - 1;`
		1731	`UNPROTECT(1);`
61696	maechler	1732	`}`
52649	ripley	1733	`}`
2	r	1734
60305	ripley	1735	`R_xlen_t iip = (R_xlen_t ) R_alloc((size_t) n, sizeof(R_xlen_t));`
41686	ripley	1736	`for (i = 0; i < n; iip[i++] = 0);`
		1737	`for (i = 0; i < n; i++)`
52649	ripley	1738	`if (pp[i] >= 0 && pp[i] < n) iip[pp[i]]++;`
		1739	`else error(_("value out of range in 'perm'"));`
41686	ripley	1740	`for (i = 0; i < n; i++)`
52649	ripley	1741	`if (iip[i] == 0) error(_("invalid '%s' argument"), "perm");`
2	r	1742
13214	maechler	1743	`/* create the stride object and permute */`
2	r	1744
60305	ripley	1745	`R_xlen_t stride = (R_xlen_t ) R_alloc((size_t) n, sizeof(R_xlen_t));`
52649	ripley	1746	`for (iip[0] = 1, i = 1; i<n; i++) iip[i] = iip[i-1] * isa[i-1];`
		1747	`for (i = 0; i < n; i++) stride[i] = iip[pp[i]];`
13214	maechler	1748
		1749	`/* also need to have the dimensions of r */`
		1750
52649	ripley	1751	`PROTECT(dimsr = allocVector(INTSXP, n));`
		1752	`int *isr = INTEGER(dimsr);`
		1753	`for (i = 0; i < n; i++) isr[i] = isa[pp[i]];`
13214	maechler	1754
		1755	`/* and away we go! iip will hold the incrementer */`
		1756
60305	ripley	1757	`R_xlen_t len = XLENGTH(a);`
1820	ihaka	1758	`PROTECT(r = allocVector(TYPEOF(a), len));`
2	r	1759
52649	ripley	1760	`for (i = 0; i < n; iip[i++] = 0);`
13214	maechler	1761
60305	ripley	1762	`R_xlen_t li, lj;`
1820	ihaka	1763	`switch (TYPEOF(a)) {`
13214	maechler	1764
1820	ihaka	1765	`case INTSXP:`
60305	ripley	1766	`for (lj = 0, li = 0; li < len; li++) {`
		1767	`INTEGER(r)[li] = INTEGER(a)[lj];`
45446	ripley	1768	`CLICKJ;`
		1769	`}`
		1770	`break;`
38546	ripley	1771
1820	ihaka	1772	`case LGLSXP:`
60305	ripley	1773	`for (lj = 0, li = 0; li < len; li++) {`
		1774	`LOGICAL(r)[li] = LOGICAL(a)[lj];`
13214	maechler	1775	`CLICKJ;`
2	r	1776	`}`
1820	ihaka	1777	`break;`
13214	maechler	1778
1820	ihaka	1779	`case REALSXP:`
60305	ripley	1780	`for (lj = 0, li = 0; li < len; li++) {`
		1781	`REAL(r)[li] = REAL(a)[lj];`
13214	maechler	1782	`CLICKJ;`
1820	ihaka	1783	`}`
		1784	`break;`
13214	maechler	1785
1820	ihaka	1786	`case CPLXSXP:`
60305	ripley	1787	`for (lj = 0, li = 0; li < len; li++) {`
		1788	`COMPLEX(r)[li].r = COMPLEX(a)[lj].r;`
		1789	`COMPLEX(r)[li].i = COMPLEX(a)[lj].i;`
13214	maechler	1790	`CLICKJ;`
1820	ihaka	1791	`}`
		1792	`break;`
13214	maechler	1793
1820	ihaka	1794	`case STRSXP:`
60305	ripley	1795	`for (lj = 0, li = 0; li < len; li++) {`
		1796	`SET_STRING_ELT(r, li, STRING_ELT(a, lj));`
13214	maechler	1797	`CLICKJ;`
1820	ihaka	1798	`}`
		1799	`break;`
13214	maechler	1800
1839	ihaka	1801	`case VECSXP:`
60305	ripley	1802	`for (lj = 0, li = 0; li < len; li++) {`
		1803	`SET_VECTOR_ELT(r, li, VECTOR_ELT(a, lj));`
13214	maechler	1804	`CLICKJ;`
1839	ihaka	1805	`}`
12976	pd	1806	`break;`
13214	maechler	1807
29902	ripley	1808	`case RAWSXP:`
60305	ripley	1809	`for (lj = 0, li = 0; li < len; li++) {`
		1810	`RAW(r)[li] = RAW(a)[lj];`
29902	ripley	1811	`CLICKJ;`
		1812	`}`
		1813	`break;`
		1814
1820	ihaka	1815	`default:`
31908	ripley	1816	`UNIMPLEMENTED_TYPE("aperm", a);`
1820	ihaka	1817	`}`
17515	ripley	1818
13214	maechler	1819	`/* handle the resize */`
52649	ripley	1820	`int resize = asLogical(CADDR(args));`
		1821	`if (resize == NA_LOGICAL) error(_("'resize' must be TRUE or FALSE"));`
13214	maechler	1822
69642	maechler	1823	`/* and handle names(dim(.)) and the dimnames if any */`
52649	ripley	1824	`if (resize) {`
69642	maechler	1825	`SEXP nmdm = getAttrib(dimsa, R_NamesSymbol);`
		1826	`if(nmdm != R_NilValue) { // dimsr needs correctly permuted names()`
		1827	`PROTECT(nmdm);`
		1828	`SEXP nm_dr = PROTECT(allocVector(STRSXP, n));`
		1829	`for (i = 0; i < n; i++) {`
		1830	`SET_STRING_ELT(nm_dr, i, STRING_ELT(nmdm, pp[i]));`
		1831	`}`
		1832	`setAttrib(dimsr, R_NamesSymbol, nm_dr);`
		1833	`UNPROTECT(2);`
		1834	`}`
		1835	`setAttrib(r, R_DimSymbol, dimsr);`
		1836
52649	ripley	1837	`PROTECT(dna = getAttrib(a, R_DimNamesSymbol));`
		1838	`if (dna != R_NilValue) {`
		1839	`SEXP dnna, dnr, dnnr;`
13214	maechler	1840
52649	ripley	1841	`PROTECT(dnr = allocVector(VECSXP, n));`
		1842	`PROTECT(dnna = getAttrib(dna, R_NamesSymbol));`
		1843	`if (dnna != R_NilValue) {`
		1844	`PROTECT(dnnr = allocVector(STRSXP, n));`
		1845	`for (i = 0; i < n; i++) {`
		1846	`SET_VECTOR_ELT(dnr, i, VECTOR_ELT(dna, pp[i]));`
		1847	`SET_STRING_ELT(dnnr, i, STRING_ELT(dnna, pp[i]));`
		1848	`}`
		1849	`setAttrib(dnr, R_NamesSymbol, dnnr);`
		1850	`UNPROTECT(1);`
		1851	`} else {`
		1852	`for (i = 0; i < n; i++)`
		1853	`SET_VECTOR_ELT(dnr, i, VECTOR_ELT(dna, pp[i]));`
		1854	`}`
		1855	`setAttrib(r, R_DimNamesSymbol, dnr);`
		1856	`UNPROTECT(2);`
13214	maechler	1857	`}`
52649	ripley	1858	`UNPROTECT(1);`
13214	maechler	1859	`}`
69642	maechler	1860	`else // !resize`
		1861	`setAttrib(r, R_DimSymbol, dimsa);`
17515	ripley	1862
52649	ripley	1863	`UNPROTECT(3); /* dimsa, r, dimsr */`
1820	ihaka	1864	`return r;`
2	r	1865	`}`
17515	ripley	1866
		1867	`/* colSums(x, n, p, na.rm) and friends */`
36990	ripley	1868	`SEXP attribute_hidden do_colsum(SEXP call, SEXP op, SEXP args, SEXP rho)`
17515	ripley	1869	`{`
		1870	`SEXP x, ans = R_NilValue;`
60282	ripley	1871	`int type;`
17521	ripley	1872	`Rboolean NaRm, keepNA;`
17515	ripley	1873
		1874	`checkArity(op, args);`
		1875	`x = CAR(args); args = CDR(args);`
64374	ripley	1876	`R_xlen_t n = asVecSize(CAR(args)); args = CDR(args);`
		1877	`R_xlen_t p = asVecSize(CAR(args)); args = CDR(args);`
17515	ripley	1878	`NaRm = asLogical(CAR(args));`
33840	ripley	1879	`if (n == NA_INTEGER \|\| n < 0)`
44512	ripley	1880	`error(_("invalid '%s' argument"), "n");`
33840	ripley	1881	`if (p == NA_INTEGER \|\| p < 0)`
44512	ripley	1882	`error(_("invalid '%s' argument"), "p");`
		1883	`if (NaRm == NA_LOGICAL) error(_("invalid '%s' argument"), "na.rm");`
17521	ripley	1884	`keepNA = !NaRm;`
17515	ripley	1885
		1886	`switch (type = TYPEOF(x)) {`
68344	maechler	1887	`case LGLSXP:`
		1888	`case INTSXP:`
17515	ripley	1889	`case REALSXP: break;`
		1890	`default:`
68344	maechler	1891	`error(_("'x' must be numeric"));`
17515	ripley	1892	`}`
68344	maechler	1893	`if (n * (double)p > XLENGTH(x))`
68923	ripley	1894	`error(_("'x' is too short")); /* PR#16367 */`
17515	ripley	1895
68344	maechler	1896	`int OP = PRIMVAL(op);`
17521	ripley	1897	`if (OP == 0 \|\| OP == 1) { /* columns */`
17515	ripley	1898	`PROTECT(ans = allocVector(REALSXP, p));`
62637	ripley	1899	`#ifdef _OPENMP`
60282	ripley	1900	`int nthreads;`
55099	ripley	1901	`/* This gives a spurious -Wunused-but-set-variable error */`
53854	luke	1902	`if (R_num_math_threads > 0)`
		1903	`nthreads = R_num_math_threads;`
		1904	`else`
		1905	`nthreads = 1; /* for now */`
		1906	`#pragma omp parallel for num_threads(nthreads) default(none) \`
59537	luke	1907	`firstprivate(x, ans, n, p, type, NaRm, keepNA, R_NaReal, R_NaInt, OP)`
53854	luke	1908	`#endif`
64374	ripley	1909	`for (R_xlen_t j = 0; j < p; j++) {`
		1910	`R_xlen_t cnt = n, i;`
60268	ripley	1911	`LDOUBLE sum = 0.0;`
17515	ripley	1912	`switch (type) {`
		1913	`case REALSXP:`
60282	ripley	1914	`{`
		1915	`double rx = REAL(x) + (R_xlen_t)nj;`
17521	ripley	1916	`if (keepNA)`
		1917	`for (sum = 0., i = 0; i < n; i++) sum += *rx++;`
		1918	`else {`
		1919	`for (cnt = 0, sum = 0., i = 0; i < n; i++, rx++)`
		1920	`if (!ISNAN(rx)) {cnt++; sum += rx;}`
60282	ripley	1921	`else if (keepNA) {sum = NA_REAL; break;} // unused`
17521	ripley	1922	`}`
17515	ripley	1923	`break;`
60282	ripley	1924	`}`
17515	ripley	1925	`case INTSXP:`
60282	ripley	1926	`{`
		1927	`int ix = INTEGER(x) + (R_xlen_t)nj;`
17521	ripley	1928	`for (cnt = 0, sum = 0., i = 0; i < n; i++, ix++)`
		1929	`if (ix != NA_INTEGER) {cnt++; sum += ix;}`
		1930	`else if (keepNA) {sum = NA_REAL; break;}`
17515	ripley	1931	`break;`
60282	ripley	1932	`}`
17515	ripley	1933	`case LGLSXP:`
60282	ripley	1934	`{`
		1935	`int ix = LOGICAL(x) + (R_xlen_t)nj;`
17521	ripley	1936	`for (cnt = 0, sum = 0., i = 0; i < n; i++, ix++)`
		1937	`if (ix != NA_LOGICAL) {cnt++; sum += ix;}`
		1938	`else if (keepNA) {sum = NA_REAL; break;}`
17515	ripley	1939	`break;`
		1940	`}`
60282	ripley	1941	`}`
59537	luke	1942	`if (OP == 1) sum /= cnt; /* gives NaN for cnt = 0 */`
59095	ripley	1943	`REAL(ans)[j] = (double) sum;`
17515	ripley	1944	`}`
		1945	`}`
59537	luke	1946	`else { /* rows */`
17515	ripley	1947	`PROTECT(ans = allocVector(REALSXP, n));`
17521	ripley	1948
59537	luke	1949	`/* allocate scratch storage to allow accumulating by columns`
		1950	`to improve cache hits */`
		1951	`int *Cnt = NULL;`
60268	ripley	1952	`LDOUBLE *rans;`
59537	luke	1953	`if(n <= 10000) {`
60268	ripley	1954	`R_CheckStack2(n * sizeof(LDOUBLE));`
		1955	`rans = (LDOUBLE ) alloca(n sizeof(LDOUBLE));`
60517	ripley	1956	`Memzero(rans, n);`
60268	ripley	1957	`} else rans = Calloc(n, LDOUBLE);`
59537	luke	1958	`if (!keepNA && OP == 3) Cnt = Calloc(n, int);`
		1959
64374	ripley	1960	`for (R_xlen_t j = 0; j < p; j++) {`
60268	ripley	1961	`LDOUBLE *ra = rans;`
59537	luke	1962	`switch (type) {`
		1963	`case REALSXP:`
60282	ripley	1964	`{`
		1965	`double rx = REAL(x) + (R_xlen_t)n j;`
17521	ripley	1966	`if (keepNA)`
64374	ripley	1967	`for (R_xlen_t i = 0; i < n; i++) ra++ += rx++;`
17521	ripley	1968	`else`
64374	ripley	1969	`for (R_xlen_t i = 0; i < n; i++, ra++, rx++)`
17521	ripley	1970	`if (!ISNAN(*rx)) {`
		1971	`ra += rx;`
59537	luke	1972	`if (OP == 3) Cnt[i]++;`
17521	ripley	1973	`}`
59537	luke	1974	`break;`
60282	ripley	1975	`}`
17515	ripley	1976	`case INTSXP:`
60282	ripley	1977	`{`
		1978	`int ix = INTEGER(x) + (R_xlen_t)n j;`
64374	ripley	1979	`for (R_xlen_t i = 0; i < n; i++, ra++, ix++)`
59537	luke	1980	`if (keepNA) {`
		1981	`if (ix != NA_INTEGER) ra += *ix;`
		1982	`else *ra = NA_REAL;`
		1983	`}`
		1984	`else if (*ix != NA_INTEGER) {`
		1985	`ra += ix;`
		1986	`if (OP == 3) Cnt[i]++;`
		1987	`}`
17515	ripley	1988	`break;`
60282	ripley	1989	`}`
17515	ripley	1990	`case LGLSXP:`
60282	ripley	1991	`{`
		1992	`int ix = LOGICAL(x) + (R_xlen_t)n j;`
64374	ripley	1993	`for (R_xlen_t i = 0; i < n; i++, ra++, ix++)`
59537	luke	1994	`if (keepNA) {`
		1995	`if (ix != NA_LOGICAL) ra += *ix;`
		1996	`else *ra = NA_REAL;`
		1997	`}`
		1998	`else if (*ix != NA_LOGICAL) {`
		1999	`ra += ix;`
		2000	`if (OP == 3) Cnt[i]++;`
		2001	`}`
17515	ripley	2002	`break;`
		2003	`}`
60282	ripley	2004	`}`
17515	ripley	2005	`}`
59537	luke	2006	`if (OP == 3) {`
		2007	`if (keepNA)`
64374	ripley	2008	`for (R_xlen_t i = 0; i < n; i++) rans[i] /= p;`
61696	maechler	2009	`else`
64374	ripley	2010	`for (R_xlen_t i = 0; i < n; i++) rans[i] /= Cnt[i];`
59537	luke	2011	`}`
64374	ripley	2012	`for (R_xlen_t i = 0; i < n; i++) REAL(ans)[i] = (double) rans[i];`
59537	luke	2013
		2014	`if (!keepNA && OP == 3) Free(Cnt);`
		2015	`if(n > 10000) Free(rans);`
17515	ripley	2016	`}`
		2017
		2018	`UNPROTECT(1);`
		2019	`return ans;`
		2020	`}`
59621	ripley	2021
		2022	`/*`
		2023	`{`
		2024	`data <- as.vector(data)`
		2025	`dim <- as.integer(dim)`
		2026	`vl <- prod(dim)`
		2027	`if (length(data) != vl) {`
68923	ripley	2028	`if (vl > .Machine$integer.max)`
		2029	`stop("'dim' specifies too large an array")`
		2030	`data <- rep(data, length.out = vl)`
59621	ripley	2031	`}`
61696	maechler	2032	`if (length(dim))`
68923	ripley	2033	`dim(data) <- dim`
61696	maechler	2034	`if (is.list(dimnames) && length(dimnames))`
68923	ripley	2035	`dimnames(data) <- dimnames`
59621	ripley	2036	`data`
		2037	`}`
		2038	`*/`
		2039
		2040	`/* array(data, dim, dimnames) */`
		2041	`SEXP attribute_hidden do_array(SEXP call, SEXP op, SEXP args, SEXP rho)`
		2042	`{`
		2043	`SEXP vals, ans, dims, dimnames;`
		2044	`R_xlen_t lendat, i, nans;`
		2045
		2046	`checkArity(op, args);`
		2047	`vals = CAR(args);`
61566	ripley	2048	`/* at least NULL can get here */`
59621	ripley	2049	`switch(TYPEOF(vals)) {`
		2050	`case LGLSXP:`
		2051	`case INTSXP:`
		2052	`case REALSXP:`
		2053	`case CPLXSXP:`
		2054	`case STRSXP:`
		2055	`case RAWSXP:`
		2056	`case EXPRSXP:`
		2057	`case VECSXP:`
		2058	`break;`
		2059	`default:`
61566	ripley	2060	`error(_("'data' must be of a vector type, was '%s'"),`
		2061	`type2char(TYPEOF(vals)));`
59621	ripley	2062	`}`
		2063	`lendat = XLENGTH(vals);`
		2064	`dims = CADR(args);`
		2065	`dimnames = CADDR(args);`
		2066	`PROTECT(dims = coerceVector(dims, INTSXP));`
		2067	`int nd = LENGTH(dims);`
59622	ripley	2068	`if (nd == 0) error(_("'dims' cannot be of length 0"));`
		2069	`double d = 1.0;`
		2070	`for (int j = 0; j < nd; j++) d *= INTEGER(dims)[j];`
59621	ripley	2071	`#ifndef LONG_VECTOR_SUPPORT`
59622	ripley	2072	`if (d > INT_MAX) error(_("too many elements specified"));`
59621	ripley	2073	`#endif`
59622	ripley	2074	`nans = (R_xlen_t) d;`
59621	ripley	2075
		2076	`PROTECT(ans = allocVector(TYPEOF(vals), nans));`
		2077	`switch(TYPEOF(vals)) {`
		2078	`case LGLSXP:`
		2079	`if (nans && lendat)`
68663	luke	2080	`xcopyLogicalWithRecycle(LOGICAL(ans), LOGICAL(vals), 0, nans,`
		2081	`lendat);`
59621	ripley	2082	`else`
		2083	`for (i = 0; i < nans; i++) LOGICAL(ans)[i] = NA_LOGICAL;`
		2084	`break;`
		2085	`case INTSXP:`
		2086	`if (nans && lendat)`
68663	luke	2087	`xcopyIntegerWithRecycle(INTEGER(ans), INTEGER(vals), 0, nans,`
		2088	`lendat);`
59621	ripley	2089	`else`
		2090	`for (i = 0; i < nans; i++) INTEGER(ans)[i] = NA_INTEGER;`
		2091	`break;`
		2092	`case REALSXP:`
		2093	`if (nans && lendat)`
68663	luke	2094	`xcopyRealWithRecycle(REAL(ans), REAL(vals), 0, nans, lendat);`
59621	ripley	2095	`else`
		2096	`for (i = 0; i < nans; i++) REAL(ans)[i] = NA_REAL;`
		2097	`break;`
		2098	`case CPLXSXP:`
		2099	`if (nans && lendat)`
68663	luke	2100	`xcopyComplexWithRecycle(COMPLEX(ans), COMPLEX(vals), 0, nans,`
		2101	`lendat);`
59621	ripley	2102	`else {`
		2103	`Rcomplex na_cmplx;`
		2104	`na_cmplx.r = NA_REAL;`
		2105	`na_cmplx.i = 0;`
		2106	`for (i = 0; i < nans; i++) COMPLEX(ans)[i] = na_cmplx;`
		2107	`}`
		2108	`break;`
		2109	`case RAWSXP:`
		2110	`if (nans && lendat)`
68663	luke	2111	`xcopyRawWithRecycle(RAW(ans), RAW(vals), 0, nans, lendat);`
59621	ripley	2112	`else`
		2113	`for (i = 0; i < nans; i++) RAW(ans)[i] = 0;`
		2114	`break;`
		2115	`case STRSXP:`
		2116	`if (nans && lendat)`
68663	luke	2117	`xcopyStringWithRecycle(ans, vals, 0, nans, lendat);`
73160	maechler	2118	`else`
		2119	`for (i = 0; i < nans; i++) SET_STRING_ELT(ans, i, NA_STRING);`
59621	ripley	2120	`break;`
73160	maechler	2121	`/* Rest are already initialized */`
59621	ripley	2122	`case VECSXP:`
		2123	`case EXPRSXP:`
66289	luke	2124	`#ifdef SWITCH_TO_REFCNT`
59621	ripley	2125	`if (nans && lendat)`
68663	luke	2126	`xcopyVectorWithRecycle(ans, vals, 0, nans, lendat);`
66289	luke	2127	`#else`
		2128	`if (nans && lendat) {`
		2129	`/* Need to guard against possible sharing of values under`
		2130	`NAMED. This is not needed with reference`
		2131	`coutning. (PR#15919) */`
		2132	`Rboolean needsmark = (lendat < nans \|\| MAYBE_REFERENCED(vals));`
		2133	`for (i = 0; i < nans; i++) {`
		2134	`SEXP elt = VECTOR_ELT(vals, i % lendat);`
		2135	`if (needsmark \|\| MAYBE_REFERENCED(elt))`
		2136	`MARK_NOT_MUTABLE(elt);`
		2137	`SET_VECTOR_ELT(ans, i, elt);`
		2138	`}`
		2139	`}`
		2140	`#endif`
59621	ripley	2141	`break;`
		2142	`default:`
		2143	`// excluded above`
		2144	`break;`
		2145	`}`
		2146
59622	ripley	2147	`ans = dimgets(ans, dims);`
69786	maechler	2148	`if(!isNull(dimnames) && length(dimnames) > 0)`
59622	ripley	2149	`ans = dimnamesgets(ans, dimnames);`
61696	maechler	2150
59621	ripley	2151	`UNPROTECT(2);`
		2152	`return ans;`
		2153	`}`
60287	ripley	2154
		2155	`SEXP attribute_hidden do_diag(SEXP call, SEXP op, SEXP args, SEXP rho)`
		2156	`{`
		2157	`SEXP ans, x, snr, snc;`
		2158	`int nr = 1, nc = 1, nprotect = 1;`
		2159
		2160	`checkArity(op, args);`
		2161	`x = CAR(args);`
		2162	`snr = CADR(args);`
		2163	`snc = CADDR(args);`
		2164	`nr = asInteger(snr);`
		2165	`if (nr == NA_INTEGER)`
		2166	`error(_("invalid 'nrow' value (too large or NA)"));`
		2167	`if (nr < 0)`
		2168	`error(_("invalid 'nrow' value (< 0)"));`
		2169	`nc = asInteger(snc);`
		2170	`if (nc == NA_INTEGER)`
		2171	`error(_("invalid 'ncol' value (too large or NA)"));`
		2172	`if (nc < 0)`
		2173	`error(_("invalid 'ncol' value (< 0)"));`
61696	maechler	2174	`int mn = (nr < nc) ? nr : nc;`
70543	murdoch	2175	`if (mn > 0 && length(x) == 0)`
60287	ripley	2176	`error(_("'x' must have positive length"));`
		2177
66495	maechler	2178	`#ifndef LONG_VECTOR_SUPPORT`
60287	ripley	2179	`if ((double)nr * (double)nc > INT_MAX)`
		2180	`error(_("too many elements specified"));`
		2181	`#endif`
		2182
70665	maechler	2183	`int nx = LENGTH(x);`
		2184	`R_xlen_t NR = nr;`
		2185
		2186	`#define mk_DIAG(_zero_) \`
		2187	`for (R_xlen_t i = 0; i < NR*nc; i++) ra[i] = _zero_; \`
		2188	`R_xlen_t i, i1; \`
		2189	`MOD_ITERATE1(mn, nx, i, i1, { \`
		2190	`ra[i * (NR+1)] = rx[i1]; \`
		2191	`});`
		2192
		2193	`switch(TYPEOF(x)) {`
		2194
		2195	`case REALSXP:`
		2196	`{`
		2197	`#define mk_REAL_DIAG \`
		2198	`PROTECT(ans = allocMatrix(REALSXP, nr, nc)); \`
		2199	`double rx = REAL(x), ra = REAL(ans); \`
		2200	`mk_DIAG(0.0)`
		2201
		2202	`mk_REAL_DIAG;`
		2203	`break;`
		2204	`}`
		2205	`case CPLXSXP:`
		2206	`{`
60287	ripley	2207	`PROTECT(ans = allocMatrix(CPLXSXP, nr, nc));`
		2208	`int nx = LENGTH(x);`
		2209	`R_xlen_t NR = nr;`
		2210	`Rcomplex rx = COMPLEX(x), ra = COMPLEX(ans), zero;`
		2211	`zero.r = zero.i = 0.0;`
70665	maechler	2212	`mk_DIAG(zero);`
		2213	`break;`
60287	ripley	2214	`}`
70665	maechler	2215	`case INTSXP:`
		2216	`{`
		2217	`PROTECT(ans = allocMatrix(INTSXP, nr, nc));`
		2218	`int rx = INTEGER(x), ra = INTEGER(ans);`
		2219	`mk_DIAG(0);`
		2220	`break;`
		2221	`}`
		2222	`case LGLSXP:`
		2223	`{`
		2224	`PROTECT(ans = allocMatrix(LGLSXP, nr, nc));`
		2225	`int rx = LOGICAL(x), ra = LOGICAL(ans);`
		2226	`mk_DIAG(0);`
		2227	`break;`
		2228	`}`
		2229	`case RAWSXP:`
		2230	`{`
		2231	`PROTECT(ans = allocMatrix(RAWSXP, nr, nc));`
		2232	`Rbyte rx = RAW(x), ra = RAW(ans);`
		2233	`mk_DIAG((Rbyte) 0);`
		2234	`break;`
		2235	`}`
		2236	`default: {`
		2237	`PROTECT(x = coerceVector(x, REALSXP));`
		2238	`nprotect++;`
		2239	`mk_REAL_DIAG;`
		2240	`}`
		2241	`}`
		2242	`#undef mk_REAL_DIAG`
		2243	`#undef mk_DIAG`
60287	ripley	2244	`UNPROTECT(nprotect);`
		2245	`return ans;`
		2246	`}`
60321	ripley	2247
		2248
60335	ripley	2249	`/* backsolve(r, b, k, upper.tri, transpose) */`
60321	ripley	2250	`SEXP attribute_hidden do_backsolve(SEXP call, SEXP op, SEXP args, SEXP rho)`
		2251	`{`
		2252	`int nprot = 1;`
		2253	`checkArity(op, args);`
		2254
		2255	`SEXP r = CAR(args); args = CDR(args);`
60335	ripley	2256	`SEXP b = CAR(args); args = CDR(args);`
		2257	`int nrr = nrows(r), nrb = nrows(b), ncb = ncols(b);`
60321	ripley	2258	`int k = asInteger(CAR(args)); args = CDR(args);`
60335	ripley	2259	`/* k is the number of rows to be used: there must be at least that`
		2260	`many rows and cols in the rhs and at least that many rows on`
		2261	`the rhs.`
		2262	`*/`
61696	maechler	2263	`if (k == NA_INTEGER \|\| k <= 0 \|\| k > nrr \|\| k > ncols(r) \|\| k > nrb)`
60321	ripley	2264	`error(_("invalid '%s' argument"), "k");`
		2265	`int upper = asLogical(CAR(args)); args = CDR(args);`
		2266	`if (upper == NA_INTEGER) error(_("invalid '%s' argument"), "upper.tri");`
		2267	`int trans = asLogical(CAR(args));`
		2268	`if (trans == NA_INTEGER) error(_("invalid '%s' argument"), "transpose");`
		2269	`if (TYPEOF(r) != REALSXP) {PROTECT(r = coerceVector(r, REALSXP)); nprot++;}`
60335	ripley	2270	`if (TYPEOF(b) != REALSXP) {PROTECT(b = coerceVector(b, REALSXP)); nprot++;}`
		2271	`double *rr = REAL(r);`
		2272
		2273	`/* check for zeros on diagonal of r: only k row/cols are used. */`
		2274	`size_t incr = nrr + 1;`
		2275	`for(int i = 0; i < k; i++) { /* check for zeros on diagonal */`
		2276	`if (rr[i * incr] == 0.0)`
		2277	`error(_("singular matrix in 'backsolve'. First zero in diagonal [%d]"),`
		2278	`i + 1);`
		2279	`}`
		2280
62618	ripley	2281	`SEXP ans = PROTECT(allocMatrix(REALSXP, k, ncb));`
60335	ripley	2282	`if (k > 0 && ncb > 0) {`
		2283	`/* copy (part) cols of b to ans */`
		2284	`for(R_xlen_t j = 0; j < ncb; j++)`
		2285	`memcpy(REAL(ans) + jk, REAL(b) + jnrb, (size_t)k *sizeof(double));`
		2286	`double one = 1.0;`
61696	maechler	2287	`F77_CALL(dtrsm)("L", upper ? "U" : "L", trans ? "T" : "N", "N",`
76585	ripley	2288	`&k, &ncb, &one, rr, &nrr, REAL(ans), &k`
		2289	`FCONE FCONE FCONE FCONE);`
60335	ripley	2290	`}`
60321	ripley	2291	`UNPROTECT(nprot);`
		2292	`return ans;`
		2293	`}`
60327	ripley	2294
		2295	`/* max.col(m, ties.method) */`
		2296	`SEXP attribute_hidden do_maxcol(SEXP call, SEXP op, SEXP args, SEXP rho)`
		2297	`{`
		2298	`checkArity(op, args);`
		2299	`SEXP m = CAR(args);`
		2300	`int method = asInteger(CADR(args));`
		2301	`int nr = nrows(m), nc = ncols(m), nprot = 1;`
62618	ripley	2302	`if (TYPEOF(m) != REALSXP) {PROTECT(m = coerceVector(m, REALSXP)); nprot++;}`
		2303	`SEXP ans = PROTECT(allocVector(INTSXP, nr));`
60327	ripley	2304	`R_max_col(REAL(m), &nr, &nc, INTEGER(ans), &method);`
		2305	`UNPROTECT(nprot);`
		2306	`return ans;`
		2307	`}`

The R Project SVN R

(root)/trunk/src/main/array.c – Rev 77458