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

The R Project SVN R

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