The R Project SVN R

Rev

Rev 62618 | Rev 64374 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 r 1 
/*
2506 maechler 2 
 *  R : A Computer Language for Statistical Data Analysis
2 r 3 
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
59035 ripley 4 
 *  Copyright (C) 1998-2012   The R Core Team
59239 ripley 5 
 *  Copyright (C) 2002-2008   The R Foundation
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
19 
 *  http://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>
35296 ripley 29 
#include <R_ext/RS.h>     /* for Calloc/Free */
13985 ripley 30 
#include <R_ext/Applic.h> /* for dgemm */
1839 ihaka 31 
 
26467 maechler 32 
/* "GetRowNames" and "GetColNames" are utility routines which
33 
 * locate and return the row names and column names from the
34 
 * dimnames attribute of a matrix.  They are useful because
35 
 * old versions of R used pair-based lists for dimnames
36 
 * whereas recent versions use vector based lists.
1881 ihaka 37 
 
26467 maechler 38 
 * These are now very old, plus
39 
 * ``When the "dimnames" attribute is
40 
 *   grabbed off an array it is always adjusted to be a vector.''
53676 ripley 41 
 
42 
 They are used in bind.c and subset.c, and advertised in Rinternals.h
26467 maechler 43 
*/
1839 ihaka 44 
SEXP GetRowNames(SEXP dimnames)
45 
{
46 
    if (TYPEOF(dimnames) == VECSXP)
10172 luke 47 
	return VECTOR_ELT(dimnames, 0);
1839 ihaka 48 
    else
1858 ihaka 49 
	return R_NilValue;
2213 maechler 50 
}
1839 ihaka 51 
 
52 
SEXP GetColNames(SEXP dimnames)
53 
{
54 
    if (TYPEOF(dimnames) == VECSXP)
10172 luke 55 
	return VECTOR_ELT(dimnames, 1);
1839 ihaka 56 
    else
2213 maechler 57 
	return R_NilValue;
1839 ihaka 58 
}
59 
 
36990 ripley 60 
SEXP attribute_hidden do_matrix(SEXP call, SEXP op, SEXP args, SEXP rho)
2 r 61 
{
53676 ripley 62 
    SEXP vals, ans, snr, snc, dimnames;
59052 ripley 63 
    int nr = 1, nc = 1, byrow, miss_nr, miss_nc;
64 
    R_xlen_t lendat;
2 r 65 
 
1820 ihaka 66 
    checkArity(op, args);
53697 ripley 67 
    vals = CAR(args); args = CDR(args);
68 
    switch(TYPEOF(vals)) {
69 
	case LGLSXP:
70 
	case INTSXP:
71 
	case REALSXP:
72 
	case CPLXSXP:
73 
	case STRSXP:
74 
	case RAWSXP:
75 
	case EXPRSXP:
76 
	case VECSXP:
77 
	    break;
78 
	default:
61566 ripley 79 
	    error(_("'data' must be of a vector type, was '%s'"),
80 
		type2char(TYPEOF(vals)));
53697 ripley 81 
    }
59052 ripley 82 
    lendat = XLENGTH(vals);
53697 ripley 83 
    snr = CAR(args); args = CDR(args);
84 
    snc = CAR(args); args = CDR(args);
85 
    byrow = asLogical(CAR(args)); args = CDR(args);
31894 ripley 86 
    if (byrow == NA_INTEGER)
44512 ripley 87 
	error(_("invalid '%s' argument"), "byrow");
53697 ripley 88 
    dimnames = CAR(args);
53900 ripley 89 
    args = CDR(args);
90 
    miss_nr = asLogical(CAR(args)); args = CDR(args);
91 
    miss_nc = asLogical(CAR(args));
2 r 92 
 
53697 ripley 93 
    if (!miss_nr) {
94 
	if (!isNumeric(snr)) error(_("non-numeric matrix extent"));
95 
	nr = asInteger(snr);
96 
	if (nr == NA_INTEGER)
97 
	    error(_("invalid 'nrow' value (too large or NA)"));
98 
	if (nr < 0)
99 
	    error(_("invalid 'nrow' value (< 0)"));
100 
    }
101 
    if (!miss_nc) {
102 
	if (!isNumeric(snc)) error(_("non-numeric matrix extent"));
103 
	nc = asInteger(snc);
104 
	if (nc == NA_INTEGER)
105 
	    error(_("invalid 'ncol' value (too large or NA)"));
106 
	if (nc < 0)
107 
	    error(_("invalid 'ncol' value (< 0)"));
108 
    }
59090 ripley 109 
    if (miss_nr && miss_nc) {
110 
	if (lendat > INT_MAX) error("data is too long");
111 
	nr = (int) lendat;
112 
    } else if (miss_nr) {
113 
	if (lendat > (double) nc * INT_MAX) error("data is too long");
59191 ripley 114 
	nr = (int) ceil((double) lendat / (double) nc);
59090 ripley 115 
    } else if (miss_nc) {
116 
	if (lendat > (double) nr * INT_MAX) error("data is too long");
59191 ripley 117 
	nc = (int) ceil((double) lendat / (double) nr);
59090 ripley 118 
    }
2 r 119 
 
59090 ripley 120 
    if(lendat > 0) {
59052 ripley 121 
	R_xlen_t nrc = (R_xlen_t) nr * nc;
122 
	if (lendat > 1 && nrc % lendat != 0) {
28014 ripley 123 
	    if (((lendat > nr) && (lendat / nr) * nr != lendat) ||
124 
		((lendat < nr) && (nr / lendat) * lendat != nr))
41686 ripley 125 
		warning(_("data length [%d] is not a sub-multiple or multiple of the number of rows [%d]"), lendat, nr);
28014 ripley 126 
	    else if (((lendat > nc) && (lendat / nc) * nc != lendat) ||
127 
		     ((lendat < nc) && (nc / lendat) * lendat != nc))
41686 ripley 128 
		warning(_("data length [%d] is not a sub-multiple or multiple of the number of columns [%d]"), lendat, nc);
28014 ripley 129 
	}
59052 ripley 130 
	else if ((lendat > 1) && (nrc == 0)){
32860 ripley 131 
	    warning(_("data length exceeds size of matrix"));
2 r 132 
	}
28014 ripley 133 
    }
2 r 134 
 
61696 maechler 135 
#ifndef LONG_VECTOR_SUPPORT
136 
    if ((double)nr * (double)nc > INT_MAX)
41680 ripley 137 
	error(_("too many elements specified"));
59052 ripley 138 
#endif
28349 ripley 139 
 
43007 ripley 140 
    PROTECT(ans = allocMatrix(TYPEOF(vals), nr, nc));
28014 ripley 141 
    if(lendat) {
142 
	if (isVector(vals))
43007 ripley 143 
	    copyMatrix(ans, vals, byrow);
28014 ripley 144 
	else
43007 ripley 145 
	    copyListMatrix(ans, vals, byrow);
28014 ripley 146 
    } else if (isVector(vals)) { /* fill with NAs */
60540 urbaneks 147 
	R_xlen_t N = (R_xlen_t) nr * nc, i;
28014 ripley 148 
	switch(TYPEOF(vals)) {
149 
	case STRSXP:
60540 urbaneks 150 
	    for (i = 0; i < N; i++)
151 
		SET_STRING_ELT(ans, i, NA_STRING);
28014 ripley 152 
	    break;
153 
	case LGLSXP:
60540 urbaneks 154 
	    for (i = 0; i < N; i++)
155 
		LOGICAL(ans)[i] = NA_LOGICAL;
28014 ripley 156 
	    break;
157 
	case INTSXP:
60540 urbaneks 158 
	    for (i = 0; i < N; i++)
159 
		INTEGER(ans)[i] = NA_INTEGER;
28014 ripley 160 
	    break;
161 
	case REALSXP:
60540 urbaneks 162 
	    for (i = 0; i < N; i++)
163 
		REAL(ans)[i] = NA_REAL;
28014 ripley 164 
	    break;
165 
	case CPLXSXP:
166 
	    {
167 
		Rcomplex na_cmplx;
168 
		na_cmplx.r = NA_REAL;
169 
		na_cmplx.i = 0;
60540 urbaneks 170 
		for (i = 0; i < N; i++)
171 
		    COMPLEX(ans)[i] = na_cmplx;
28014 ripley 172 
	    }
173 
	    break;
29902 ripley 174 
	case RAWSXP:
60540 urbaneks 175 
	    memset(RAW(ans), 0, N);
29902 ripley 176 
	    break;
32350 maechler 177 
	default:
31875 ripley 178 
	    /* don't fill with anything */
179 
	    ;
28014 ripley 180 
	}
181 
    }
48092 maechler 182 
    if(!isNull(dimnames)&& length(dimnames) > 0)
47049 ripley 183 
	ans = dimnamesgets(ans, dimnames);
1820 ihaka 184 
    UNPROTECT(1);
43007 ripley 185 
    return ans;
2 r 186 
}
187 
 
188 
 
189 
SEXP allocMatrix(SEXPTYPE mode, int nrow, int ncol)
190 
{
1820 ihaka 191 
    SEXP s, t;
59052 ripley 192 
    R_xlen_t n;
2 r 193 
 
1820 ihaka 194 
    if (nrow < 0 || ncol < 0)
32860 ripley 195 
	error(_("negative extents to matrix"));
59052 ripley 196 
#ifndef LONG_VECTOR_SUPPORT
28349 ripley 197 
    if ((double)nrow * (double)ncol > INT_MAX)
32860 ripley 198 
	error(_("allocMatrix: too many elements specified"));
59052 ripley 199 
#endif
200 
    n = ((R_xlen_t) nrow) * ncol;
1820 ihaka 201 
    PROTECT(s = allocVector(mode, n));
202 
    PROTECT(t = allocVector(INTSXP, 2));
203 
    INTEGER(t)[0] = nrow;
204 
    INTEGER(t)[1] = ncol;
205 
    setAttrib(s, R_DimSymbol, t);
206 
    UNPROTECT(2);
207 
    return s;
2 r 208 
}
209 
 
42325 bates 210 
/**
211 
 * Allocate a 3-dimensional array
212 
 *
213 
 * @param mode The R mode (e.g. INTSXP)
214 
 * @param nrow number of rows
215 
 * @param ncol number of columns
216 
 * @param nface number of faces
217 
 *
218 
 * @return A 3-dimensional array of the indicated dimensions and mode
219 
 */
220 
SEXP alloc3DArray(SEXPTYPE mode, int nrow, int ncol, int nface)
221 
{
222 
    SEXP s, t;
59052 ripley 223 
    R_xlen_t n;
1130 maechler 224 
 
42325 bates 225 
    if (nrow < 0 || ncol < 0 || nface < 0)
226 
	error(_("negative extents to 3D array"));
59052 ripley 227 
#ifndef LONG_VECTOR_SUPPORT
42325 bates 228 
    if ((double)nrow * (double)ncol * (double)nface > INT_MAX)
60844 ripley 229 
	error(_("'alloc3Darray': too many elements specified"));
59052 ripley 230 
#endif
231 
    n = ((R_xlen_t) nrow) * ncol * nface;
42325 bates 232 
    PROTECT(s = allocVector(mode, n));
233 
    PROTECT(t = allocVector(INTSXP, 3));
234 
    INTEGER(t)[0] = nrow;
235 
    INTEGER(t)[1] = ncol;
236 
    INTEGER(t)[2] = nface;
237 
    setAttrib(s, R_DimSymbol, t);
238 
    UNPROTECT(2);
239 
    return s;
240 
}
241 
 
242 
 
2 r 243 
SEXP allocArray(SEXPTYPE mode, SEXP dims)
244 
{
1820 ihaka 245 
    SEXP array;
59052 ripley 246 
    int i;
59090 ripley 247 
    R_xlen_t n = 1;
248 
    double dn = 1;
2 r 249 
 
28349 ripley 250 
    for (i = 0; i < LENGTH(dims); i++) {
251 
	dn *= INTEGER(dims)[i];
59052 ripley 252 
#ifndef LONG_VECTOR_SUPPORT
28349 ripley 253 
	if(dn > INT_MAX)
60844 ripley 254 
	    error(_("'allocArray': too many elements specified by 'dims'"));
59052 ripley 255 
#endif
28349 ripley 256 
	n *= INTEGER(dims)[i];
257 
    }
2 r 258 
 
1820 ihaka 259 
    PROTECT(dims = duplicate(dims));
260 
    PROTECT(array = allocVector(mode, n));
261 
    setAttrib(array, R_DimSymbol, dims);
262 
    UNPROTECT(2);
263 
    return array;
2 r 264 
}
265 
 
1820 ihaka 266 
/* DropDims strips away redundant dimensioning information. */
267 
/* If there is an appropriate dimnames attribute the correct */
268 
/* element is extracted and attached to the vector as a names */
269 
/* attribute.  Note that this function mutates x. */
270 
/* Duplication should occur before this is called. */
2 r 271 
 
272 
SEXP DropDims(SEXP x)
273 
{
43878 ripley 274 
    SEXP dims, dimnames, newnames = R_NilValue;
1839 ihaka 275 
    int i, n, ndims;
2 r 276 
 
1820 ihaka 277 
    PROTECT(x);
278 
    dims = getAttrib(x, R_DimSymbol);
279 
    dimnames = getAttrib(x, R_DimNamesSymbol);
2 r 280 
 
1820 ihaka 281 
    /* Check that dropping will actually do something. */
282 
    /* (1) Check that there is a "dim" attribute. */
2 r 283 
 
1839 ihaka 284 
    if (dims == R_NilValue) {
1820 ihaka 285 
	UNPROTECT(1);
286 
	return x;
287 
    }
1839 ihaka 288 
    ndims = LENGTH(dims);
2 r 289 
 
1839 ihaka 290 
    /* (2) Check whether there are redundant extents */
1820 ihaka 291 
    n = 0;
1839 ihaka 292 
    for (i = 0; i < ndims; i++)
293 
	if (INTEGER(dims)[i] != 1) n++;
294 
    if (n == ndims) {
1820 ihaka 295 
	UNPROTECT(1);
296 
	return x;
297 
    }
2 r 298 
 
1839 ihaka 299 
    if (n <= 1) {
43878 ripley 300 
	/* We have reduced to a vector result.
301 
	   If that has length one, it is ambiguous which dimnames to use,
302 
	   so use it if there is only one (as from R 2.7.0).
303 
	 */
1820 ihaka 304 
	if (dimnames != R_NilValue) {
59052 ripley 305 
	    if(XLENGTH(x) != 1) {
43878 ripley 306 
		for (i = 0; i < LENGTH(dims); i++) {
1839 ihaka 307 
		    if (INTEGER(dims)[i] != 1) {
10172 luke 308 
			newnames = VECTOR_ELT(dimnames, i);
1839 ihaka 309 
			break;
310 
		    }
2 r 311 
		}
43878 ripley 312 
	    } else { /* drop all dims: keep names if unambiguous */
313 
		int cnt;
314 
		for(i = 0, cnt = 0; i < LENGTH(dims); i++)
315 
		    if(VECTOR_ELT(dimnames, i) != R_NilValue) cnt++;
316 
		if(cnt == 1)
317 
		    for (i = 0; i < LENGTH(dims); i++) {
318 
			newnames = VECTOR_ELT(dimnames, i);
319 
			if(newnames != R_NilValue) break;
1839 ihaka 320 
		    }
1820 ihaka 321 
	    }
2 r 322 
	}
1820 ihaka 323 
	PROTECT(newnames);
324 
	setAttrib(x, R_DimNamesSymbol, R_NilValue);
325 
	setAttrib(x, R_DimSymbol, R_NilValue);
326 
	setAttrib(x, R_NamesSymbol, newnames);
48062 jmc 327 
	/* FIXME: the following is desirable, but pointless as long as
328 
	   subset.c & others have a contrary version that leaves the
329 
	   S4 class in, incorrectly, in the case of vectors.  JMC
330 
	   3/3/09 */
331 
/* 	if(IS_S4_OBJECT(x)) {/\* no longer valid subclass of array or
332 
 	matrix *\/ */
333 
/* 	    setAttrib(x, R_ClassSymbol, R_NilValue); */
334 
/* 	    UNSET_S4_OBJECT(x); */
335 
/* 	} */
1820 ihaka 336 
	UNPROTECT(1);
43878 ripley 337 
    } else {
1839 ihaka 338 
	/* We have a lower dimensional array. */
6994 pd 339 
	SEXP newdims, dnn, newnamesnames = R_NilValue;
340 
	dnn = getAttrib(dimnames, R_NamesSymbol);
1820 ihaka 341 
	PROTECT(newdims = allocVector(INTSXP, n));
1858 ihaka 342 
	for (i = 0, n = 0; i < ndims; i++)
1839 ihaka 343 
	    if (INTEGER(dims)[i] != 1)
1820 ihaka 344 
		INTEGER(newdims)[n++] = INTEGER(dims)[i];
6994 pd 345 
	if (!isNull(dimnames)) {
1858 ihaka 346 
	    int havenames = 0;
347 
	    for (i = 0; i < ndims; i++)
10172 luke 348 
		if (INTEGER(dims)[i] != 1 &&
349 
		    VECTOR_ELT(dimnames, i) != R_NilValue)
1858 ihaka 350 
		    havenames = 1;
351 
	    if (havenames) {
2715 pd 352 
		PROTECT(newnames = allocVector(VECSXP, n));
6994 pd 353 
		PROTECT(newnamesnames = allocVector(STRSXP, n));
354 
		for (i = 0, n = 0; i < ndims; i++) {
355 
		    if (INTEGER(dims)[i] != 1) {
356 
			if(!isNull(dnn))
10172 luke 357 
			    SET_STRING_ELT(newnamesnames, n,
358 
					   STRING_ELT(dnn, i));
359 
			SET_VECTOR_ELT(newnames, n++, VECTOR_ELT(dimnames, i));
6994 pd 360 
		    }
1858 ihaka 361 
		}
1820 ihaka 362 
	    }
1858 ihaka 363 
	    else dimnames = R_NilValue;
1839 ihaka 364 
	}
1858 ihaka 365 
	PROTECT(dimnames);
1820 ihaka 366 
	setAttrib(x, R_DimNamesSymbol, R_NilValue);
367 
	setAttrib(x, R_DimSymbol, newdims);
2506 maechler 368 
	if (dimnames != R_NilValue)
2715 pd 369 
	{
6994 pd 370 
	    if(!isNull(dnn))
371 
		setAttrib(newnames, R_NamesSymbol, newnamesnames);
2506 maechler 372 
	    setAttrib(x, R_DimNamesSymbol, newnames);
6994 pd 373 
	    UNPROTECT(2);
2715 pd 374 
	}
1858 ihaka 375 
	UNPROTECT(2);
1820 ihaka 376 
    }
377 
    UNPROTECT(1);
378 
    return x;
2 r 379 
}
380 
 
36990 ripley 381 
SEXP attribute_hidden do_drop(SEXP call, SEXP op, SEXP args, SEXP rho)
2 r 382 
{
1820 ihaka 383 
    SEXP x, xdims;
384 
    int i, n, shorten;
2 r 385 
 
1820 ihaka 386 
    checkArity(op, args);
387 
    x = CAR(args);
1839 ihaka 388 
    if ((xdims = getAttrib(x, R_DimSymbol)) != R_NilValue) {
1820 ihaka 389 
	n = LENGTH(xdims);
390 
	shorten = 0;
1839 ihaka 391 
	for (i = 0; i < n; i++)
392 
	    if (INTEGER(xdims)[i] == 1) shorten = 1;
393 
	if (shorten) {
394 
	    if (NAMED(x)) x = duplicate(x);
1820 ihaka 395 
	    x = DropDims(x);
2 r 396 
	}
1820 ihaka 397 
    }
398 
    return x;
2 r 399 
}
400 
 
1820 ihaka 401 
/* Length of Primitive Objects */
2 r 402 
 
36990 ripley 403 
SEXP attribute_hidden do_length(SEXP call, SEXP op, SEXP args, SEXP rho)
2 r 404 
{
40176 ripley 405 
    checkArity(op, args);
51267 ripley 406 
    check1arg(args, call, "x");
15783 rgentlem 407 
 
60156 ripley 408 
    SEXP x = CAR(args), ans;
59185 ripley 409 
 
60156 ripley 410 
    if (isObject(x) &&
59239 ripley 411 
       DispatchOrEval(call, op, "length", args, rho, &ans, 0, 1)) {
412 
	if (length(ans) == 1 && TYPEOF(ans) == REALSXP) {
59241 ripley 413 
	    double d = REAL(ans)[0];
59249 ripley 414 
	    if (R_FINITE(d) && d >= 0. && d <= INT_MAX && floor(d) == d)
59241 ripley 415 
		return coerceVector(ans, INTSXP);
59239 ripley 416 
	}
59224 ripley 417 
	return(ans);
59239 ripley 418 
    }
17515 ripley 419 
 
59185 ripley 420 
#ifdef LONG_VECTOR_SUPPORT
60132 ripley 421 
    // or use IS_LONG_VEC
59185 ripley 422 
    R_xlen_t len = xlength(x);
60132 ripley 423 
    if (len > INT_MAX) return ScalarReal((double) len);
424 
#endif
59185 ripley 425 
    return ScalarInteger(length(x));
2 r 426 
}
427 
 
428 
 
36990 ripley 429 
SEXP attribute_hidden do_rowscols(SEXP call, SEXP op, SEXP args, SEXP rho)
2 r 430 
{
41762 ripley 431 
    SEXP x, ans;
1820 ihaka 432 
    int i, j, nr, nc;
2 r 433 
 
41712 ripley 434 
    checkArity(op, args);
59052 ripley 435 
    /* This is the dimensions vector */
41762 ripley 436 
    x = CAR(args);
437 
    if (!isInteger(x) || LENGTH(x) != 2)
60844 ripley 438 
	error(_("a matrix-like object is required as argument to '%s'"),
439 
	      (PRIMVAL(op) == 2) ? "col" : "row");
2 r 440 
 
41762 ripley 441 
    nr = INTEGER(x)[0];
442 
    nc = INTEGER(x)[1];
2 r 443 
 
1820 ihaka 444 
    ans = allocMatrix(INTSXP, nr, nc);
2 r 445 
 
59052 ripley 446 
    R_xlen_t NR = nr;
1820 ihaka 447 
    switch (PRIMVAL(op)) {
448 
    case 1:
449 
	for (i = 0; i < nr; i++)
450 
	    for (j = 0; j < nc; j++)
59052 ripley 451 
		INTEGER(ans)[i + j * NR] = i + 1;
1820 ihaka 452 
	break;
453 
    case 2:
454 
	for (i = 0; i < nr; i++)
455 
	    for (j = 0; j < nc; j++)
59052 ripley 456 
		INTEGER(ans)[i + j * NR] = j + 1;
1820 ihaka 457 
	break;
458 
    }
459 
    return ans;
2 r 460 
}
461 
 
13942 bates 462 
static void matprod(double *x, int nrx, int ncx,
463 
		    double *y, int nry, int ncy, double *z)
464 
{
465 
    char *transa = "N", *transb = "N";
37393 ripley 466 
    double one = 1.0, zero = 0.0;
60268 ripley 467 
    LDOUBLE sum;
27297 ripley 468 
    Rboolean have_na = FALSE;
60241 ripley 469 
    R_xlen_t NRX = nrx, NRY = nry;
27297 ripley 470 
 
13996 duncan 471 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
27297 ripley 472 
	/* Don't trust the BLAS to handle NA/NaNs correctly: PR#4582
59068 ripley 473 
	 * The test is only O(n) here.
27297 ripley 474 
	 */
59068 ripley 475 
	for (R_xlen_t i = 0; i < NRX*ncx; i++)
27297 ripley 476 
	    if (ISNAN(x[i])) {have_na = TRUE; break;}
32350 maechler 477 
	if (!have_na)
59068 ripley 478 
	    for (R_xlen_t i = 0; i < NRY*ncy; i++)
27297 ripley 479 
		if (ISNAN(y[i])) {have_na = TRUE; break;}
480 
	if (have_na) {
60241 ripley 481 
	    for (int i = 0; i < nrx; i++)
482 
		for (int k = 0; k < ncy; k++) {
27297 ripley 483 
		    sum = 0.0;
60241 ripley 484 
		    for (int j = 0; j < ncx; j++)
59068 ripley 485 
			sum += x[i + j * NRX] * y[j + k * NRY];
60241 ripley 486 
		    z[i + k * NRX] = (double) sum;
27297 ripley 487 
		}
488 
	} else
489 
	    F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
490 
			    x, &nrx, y, &nry, &zero, z, &nrx);
491 
    } else /* zero-extent operations should return zeroes */
60241 ripley 492 
	for(R_xlen_t i = 0; i < NRX*ncy; i++) z[i] = 0;
27297 ripley 493 
}
571 ihaka 494 
 
6994 pd 495 
static void cmatprod(Rcomplex *x, int nrx, int ncx,
23175 ripley 496 
		     Rcomplex *y, int nry, int ncy, Rcomplex *z)
2 r 497 
{
35212 ripley 498 
#ifdef HAVE_FORTRAN_DOUBLE_COMPLEX
23175 ripley 499 
    char *transa = "N", *transb = "N";
500 
    Rcomplex one, zero;
501 
 
502 
    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
503 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
45446 ripley 504 
	F77_CALL(zgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
23175 ripley 505 
			x, &nrx, y, &nry, &zero, z, &nrx);
506 
    } else { /* zero-extent operations should return zeroes */
60241 ripley 507 
	R_xlen_t NRX = nrx;
508 
	for(R_xlen_t i = 0; i < NRX*ncy; i++) z[i].r = z[i].i = 0;
23175 ripley 509 
    }
510 
#else
1820 ihaka 511 
    int i, j, k;
37393 ripley 512 
    double xij_r, xij_i, yjk_r, yjk_i;
60268 ripley 513 
    LDOUBLE sum_i, sum_r;
2 r 514 
 
59068 ripley 515 
    R_xlen_t NRX = nrx, NRY = nry;
1839 ihaka 516 
    for (i = 0; i < nrx; i++)
517 
	for (k = 0; k < ncy; k++) {
59068 ripley 518 
	    z[i + k * NRX].r = NA_REAL;
519 
	    z[i + k * NRX].i = NA_REAL;
1820 ihaka 520 
	    sum_r = 0.0;
521 
	    sum_i = 0.0;
1839 ihaka 522 
	    for (j = 0; j < ncx; j++) {
59068 ripley 523 
		xij_r = x[i + j * NRX].r;
524 
		xij_i = x[i + j * NRX].i;
525 
		yjk_r = y[j + k * NRY].r;
526 
		yjk_i = y[j + k * NRY].i;
1820 ihaka 527 
		if (ISNAN(xij_r) || ISNAN(xij_i)
528 
		    || ISNAN(yjk_r) || ISNAN(yjk_i))
529 
		    goto next_ik;
530 
		sum_r += (xij_r * yjk_r - xij_i * yjk_i);
531 
		sum_i += (xij_r * yjk_i + xij_i * yjk_r);
532 
	    }
59068 ripley 533 
	    z[i + k * NRX].r = sum_r;
534 
	    z[i + k * NRX].i = sum_i;
1820 ihaka 535 
	next_ik:
536 
	    ;
23175 ripley 537 
	}
1016 maechler 538 
#endif
2 r 539 
}
540 
 
18753 ripley 541 
static void symcrossprod(double *x, int nr, int nc, double *z)
542 
{
543 
    char *trans = "T", *uplo = "U";
544 
    double one = 1.0, zero = 0.0;
60241 ripley 545 
    R_xlen_t NC = nc;
18753 ripley 546 
    if (nr > 0 && nc > 0) {
45446 ripley 547 
	F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc);
60241 ripley 548 
	for (int i = 1; i < nc; i++)
549 
	    for (int j = 0; j < i; j++) z[i + NC *j] = z[j + NC * i];
26666 ripley 550 
    } else { /* zero-extent operations should return zeroes */
60241 ripley 551 
	for(R_xlen_t i = 0; i < NC*NC; i++) z[i] = 0;
18753 ripley 552 
    }
26666 ripley 553 
 
18753 ripley 554 
}
555 
 
1820 ihaka 556 
static void crossprod(double *x, int nrx, int ncx,
557 
		      double *y, int nry, int ncy, double *z)
2 r 558 
{
13942 bates 559 
    char *transa = "T", *transb = "N";
560 
    double one = 1.0, zero = 0.0;
13996 duncan 561 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
45446 ripley 562 
	F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
23175 ripley 563 
			x, &nrx, y, &nry, &zero, z, &ncx);
36679 ripley 564 
    } else { /* zero-extent operations should return zeroes */
60241 ripley 565 
	R_xlen_t NCX = ncx;
566 
	for(R_xlen_t i = 0; i < NCX*ncy; i++) z[i] = 0;
26666 ripley 567 
    }
2 r 568 
}
569 
 
6994 pd 570 
static void ccrossprod(Rcomplex *x, int nrx, int ncx,
571 
		       Rcomplex *y, int nry, int ncy, Rcomplex *z)
2 r 572 
{
23175 ripley 573 
    char *transa = "T", *transb = "N";
574 
    Rcomplex one, zero;
575 
 
576 
    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
577 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
45446 ripley 578 
	F77_CALL(zgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
23175 ripley 579 
			x, &nrx, y, &nry, &zero, z, &ncx);
36679 ripley 580 
    } else { /* zero-extent operations should return zeroes */
60241 ripley 581 
	R_xlen_t NCX = ncx;
582 
	for(R_xlen_t i = 0; i < NCX*ncy; i++) z[i].r = z[i].i = 0;
26666 ripley 583 
    }
2 r 584 
}
36679 ripley 585 
 
586 
static void symtcrossprod(double *x, int nr, int nc, double *z)
587 
{
588 
    char *trans = "N", *uplo = "U";
589 
    double one = 1.0, zero = 0.0;
590 
    if (nr > 0 && nc > 0) {
45446 ripley 591 
	F77_CALL(dsyrk)(uplo, trans, &nr, &nc, &one, x, &nr, &zero, z, &nr);
60241 ripley 592 
	for (int i = 1; i < nr; i++)
593 
	    for (int j = 0; j < i; j++) z[i + nr *j] = z[j + nr * i];
36679 ripley 594 
    } else { /* zero-extent operations should return zeroes */
60241 ripley 595 
	R_xlen_t NR = nr;
596 
	for(R_xlen_t i = 0; i < NR*NR; i++) z[i] = 0;
36679 ripley 597 
    }
598 
 
599 
}
600 
 
601 
static void tcrossprod(double *x, int nrx, int ncx,
602 
		      double *y, int nry, int ncy, double *z)
603 
{
604 
    char *transa = "N", *transb = "T";
605 
    double one = 1.0, zero = 0.0;
606 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
45446 ripley 607 
	F77_CALL(dgemm)(transa, transb, &nrx, &nry, &ncx, &one,
36679 ripley 608 
			x, &nrx, y, &nry, &zero, z, &nrx);
609 
    } else { /* zero-extent operations should return zeroes */
60241 ripley 610 
	R_xlen_t NRX = nrx;
611 
	for(R_xlen_t i = 0; i < NRX*nry; i++) z[i] = 0;
36679 ripley 612 
    }
613 
}
614 
 
615 
static void tccrossprod(Rcomplex *x, int nrx, int ncx,
616 
			Rcomplex *y, int nry, int ncy, Rcomplex *z)
617 
{
618 
    char *transa = "N", *transb = "T";
619 
    Rcomplex one, zero;
620 
 
621 
    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
622 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
45446 ripley 623 
	F77_CALL(zgemm)(transa, transb, &nrx, &nry, &ncx, &one,
36679 ripley 624 
			x, &nrx, y, &nry, &zero, z, &nrx);
625 
    } else { /* zero-extent operations should return zeroes */
60241 ripley 626 
	R_xlen_t NRX = nrx;
627 
	for(R_xlen_t i = 0; i < NRX*nry; i++) z[i].r = z[i].i = 0;
36679 ripley 628 
    }
629 
}
630 
 
631 
 
632 
/* "%*%" (op = 0), crossprod (op = 1) or tcrossprod (op = 2) */
36990 ripley 633 
SEXP attribute_hidden do_matprod(SEXP call, SEXP op, SEXP args, SEXP rho)
2 r 634 
{
1820 ihaka 635 
    int ldx, ldy, nrx, ncx, nry, ncy, mode;
18753 ripley 636 
    SEXP x = CAR(args), y = CADR(args), xdims, ydims, ans;
637 
    Rboolean sym;
2 r 638 
 
59239 ripley 639 
    if (PRIMVAL(op) == 0 && /* %*% is primitive, the others are .Internal() */
45446 ripley 640 
       (IS_S4_OBJECT(x) || IS_S4_OBJECT(y))
41281 ripley 641 
       && R_has_methods(op)) {
41365 ripley 642 
	SEXP s, value;
643 
	/* Remove argument names to ensure positional matching */
644 
	for(s = args; s != R_NilValue; s = CDR(s)) SET_TAG(s, R_NilValue);
42292 ripley 645 
	value = R_possible_dispatch(call, op, args, rho, FALSE);
59239 ripley 646 
	if (value) return value;
21875 jmc 647 
    }
648 
 
18753 ripley 649 
    sym = isNull(y);
36679 ripley 650 
    if (sym && (PRIMVAL(op) > 0)) y = x;
18753 ripley 651 
    if ( !(isNumeric(x) || isComplex(x)) || !(isNumeric(y) || isComplex(y)) )
48092 maechler 652 
	errorcall(call, _("requires numeric/complex matrix/vector arguments"));
2 r 653 
 
1820 ihaka 654 
    xdims = getAttrib(x, R_DimSymbol);
655 
    ydims = getAttrib(y, R_DimSymbol);
656 
    ldx = length(xdims);
657 
    ldy = length(ydims);
2 r 658 
 
26472 ihaka 659 
    if (ldx != 2 && ldy != 2) {		/* x and y non-matrices */
1839 ihaka 660 
	if (PRIMVAL(op) == 0) {
1820 ihaka 661 
	    nrx = 1;
662 
	    ncx = LENGTH(x);
2 r 663 
	}
1820 ihaka 664 
	else {
665 
	    nrx = LENGTH(x);
666 
	    ncx = 1;
2 r 667 
	}
1820 ihaka 668 
	nry = LENGTH(y);
669 
	ncy = 1;
670 
    }
26472 ihaka 671 
    else if (ldx != 2) {		/* x not a matrix */
1820 ihaka 672 
	nry = INTEGER(ydims)[0];
673 
	ncy = INTEGER(ydims)[1];
674 
	nrx = 0;
675 
	ncx = 0;
1839 ihaka 676 
	if (PRIMVAL(op) == 0) {
26472 ihaka 677 
	    if (LENGTH(x) == nry) {	/* x as row vector */
1820 ihaka 678 
		nrx = 1;
48092 maechler 679 
		ncx = nry; /* == LENGTH(x) */
2213 maechler 680 
	    }
26472 ihaka 681 
	    else if (nry == 1) {	/* x as col vector */
1820 ihaka 682 
		nrx = LENGTH(x);
683 
		ncx = 1;
684 
	    }
2 r 685 
	}
48574 maechler 686 
	else if (PRIMVAL(op) == 1) { /* crossprod() */
36679 ripley 687 
	    if (LENGTH(x) == nry) {	/* x is a col vector */
48092 maechler 688 
		nrx = nry; /* == LENGTH(x) */
1820 ihaka 689 
		ncx = 1;
690 
	    }
50989 maechler 691 
	    /* else if (nry == 1) ... not being too tolerant
692 
	       to treat x as row vector, as t(x) *is* row vector */
2 r 693 
	}
48574 maechler 694 
	else { /* tcrossprod */
695 
	    if (LENGTH(x) == ncy) {	/* x as row vector */
696 
		nrx = 1;
697 
		ncx = ncy; /* == LENGTH(x) */
698 
	    }
50989 maechler 699 
	    else if (ncy == 1) {	/* x as col vector */
700 
		nrx = LENGTH(x);
701 
		ncx = 1;
702 
	    }
48574 maechler 703 
	}
1820 ihaka 704 
    }
26472 ihaka 705 
    else if (ldy != 2) {		/* y not a matrix */
1820 ihaka 706 
	nrx = INTEGER(xdims)[0];
707 
	ncx = INTEGER(xdims)[1];
708 
	nry = 0;
709 
	ncy = 0;
1839 ihaka 710 
	if (PRIMVAL(op) == 0) {
26472 ihaka 711 
	    if (LENGTH(y) == ncx) {	/* y as col vector */
50989 maechler 712 
		nry = ncx;
1820 ihaka 713 
		ncy = 1;
714 
	    }
26472 ihaka 715 
	    else if (ncx == 1) {	/* y as row vector */
1820 ihaka 716 
		nry = 1;
717 
		ncy = LENGTH(y);
718 
	    }
2 r 719 
	}
50989 maechler 720 
	else if (PRIMVAL(op) == 1) { /* crossprod() */
36679 ripley 721 
	    if (LENGTH(y) == nrx) {	/* y is a col vector */
50989 maechler 722 
		nry = nrx;
1820 ihaka 723 
		ncy = 1;
724 
	    }
2 r 725 
	}
50989 maechler 726 
	else { /* tcrossprod --		y is a col vector */
727 
	    nry = LENGTH(y);
728 
	    ncy = 1;
729 
	}
1820 ihaka 730 
    }
26472 ihaka 731 
    else {				/* x and y matrices */
1820 ihaka 732 
	nrx = INTEGER(xdims)[0];
733 
	ncx = INTEGER(xdims)[1];
734 
	nry = INTEGER(ydims)[0];
735 
	ncy = INTEGER(ydims)[1];
736 
    }
32350 maechler 737 
    /* nr[ow](.) and nc[ol](.) are now defined for x and y */
2 r 738 
 
1839 ihaka 739 
    if (PRIMVAL(op) == 0) {
41686 ripley 740 
	/* primitive, so use call */
1839 ihaka 741 
	if (ncx != nry)
32860 ripley 742 
	    errorcall(call, _("non-conformable arguments"));
1820 ihaka 743 
    }
36679 ripley 744 
    else if (PRIMVAL(op) == 1) {
1839 ihaka 745 
	if (nrx != nry)
41686 ripley 746 
	    error(_("non-conformable arguments"));
1820 ihaka 747 
    }
36679 ripley 748 
    else {
749 
	if (ncx != ncy)
41686 ripley 750 
	    error(_("non-conformable arguments"));
36679 ripley 751 
    }
1820 ihaka 752 
 
1839 ihaka 753 
    if (isComplex(CAR(args)) || isComplex(CADR(args)))
1820 ihaka 754 
	mode = CPLXSXP;
755 
    else
756 
	mode = REALSXP;
10172 luke 757 
    SETCAR(args, coerceVector(CAR(args), mode));
758 
    SETCADR(args, coerceVector(CADR(args), mode));
1820 ihaka 759 
 
32350 maechler 760 
    if (PRIMVAL(op) == 0) {			/* op == 0 : matprod() */
26472 ihaka 761 
 
1820 ihaka 762 
	PROTECT(ans = allocMatrix(mode, nrx, ncy));
1839 ihaka 763 
	if (mode == CPLXSXP)
1820 ihaka 764 
	    cmatprod(COMPLEX(CAR(args)), nrx, ncx,
765 
		     COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));
2 r 766 
	else
1820 ihaka 767 
	    matprod(REAL(CAR(args)), nrx, ncx,
768 
		    REAL(CADR(args)), nry, ncy, REAL(ans));
26472 ihaka 769 
 
1820 ihaka 770 
	PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));
771 
	PROTECT(ydims = getAttrib(CADR(args), R_DimNamesSymbol));
26472 ihaka 772 
 
1820 ihaka 773 
	if (xdims != R_NilValue || ydims != R_NilValue) {
32350 maechler 774 
	    SEXP dimnames, dimnamesnames, dnx=R_NilValue, dny=R_NilValue;
775 
 
776 
	    /* allocate dimnames and dimnamesnames */
777 
 
6994 pd 778 
	    PROTECT(dimnames = allocVector(VECSXP, 2));
779 
	    PROTECT(dimnamesnames = allocVector(STRSXP, 2));
780 
	    if (xdims != R_NilValue) {
32350 maechler 781 
		if (ldx == 2 || ncx == 1) {
10172 luke 782 
		    SET_VECTOR_ELT(dimnames, 0, VECTOR_ELT(xdims, 0));
32350 maechler 783 
		    dnx = getAttrib(xdims, R_NamesSymbol);
784 
		    if(!isNull(dnx))
785 
			SET_STRING_ELT(dimnamesnames, 0, STRING_ELT(dnx, 0));
8628 pd 786 
		}
6994 pd 787 
	    }
32350 maechler 788 
 
32356 maechler 789 
#define YDIMS_ET_CETERA							\
790 
	    if (ydims != R_NilValue) {					\
791 
		if (ldy == 2) {						\
792 
		    SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 1));	\
793 
		    dny = getAttrib(ydims, R_NamesSymbol);		\
794 
		    if(!isNull(dny))					\
795 
			SET_STRING_ELT(dimnamesnames, 1, STRING_ELT(dny, 1)); \
796 
		} else if (nry == 1) {					\
797 
		    SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 0));	\
798 
		    dny = getAttrib(ydims, R_NamesSymbol);		\
799 
		    if(!isNull(dny))					\
800 
			SET_STRING_ELT(dimnamesnames, 1, STRING_ELT(dny, 0)); \
801 
		}							\
802 
	    }								\
803 
									\
804 
	    /* We sometimes attach a dimnames attribute			\
805 
	     * whose elements are all NULL ...				\
806 
	     * This is ugly but causes no real damage.			\
807 
	     * Now (2.1.0 ff), we don't anymore: */			\
808 
	    if (VECTOR_ELT(dimnames,0) != R_NilValue ||			\
809 
		VECTOR_ELT(dimnames,1) != R_NilValue) {			\
810 
		if (dnx != R_NilValue || dny != R_NilValue)		\
811 
		    setAttrib(dimnames, R_NamesSymbol, dimnamesnames);	\
812 
		setAttrib(ans, R_DimNamesSymbol, dimnames);		\
813 
	    }								\
32350 maechler 814 
	    UNPROTECT(2)
815 
 
816 
	    YDIMS_ET_CETERA;
2 r 817 
	}
1820 ihaka 818 
    }
26472 ihaka 819 
 
36679 ripley 820 
    else if (PRIMVAL(op) == 1) {	/* op == 1: crossprod() */
26472 ihaka 821 
 
1820 ihaka 822 
	PROTECT(ans = allocMatrix(mode, ncx, ncy));
1839 ihaka 823 
	if (mode == CPLXSXP)
23175 ripley 824 
	    if(sym)
825 
		ccrossprod(COMPLEX(CAR(args)), nrx, ncx,
826 
			   COMPLEX(CAR(args)), nry, ncy, COMPLEX(ans));
827 
	    else
828 
		ccrossprod(COMPLEX(CAR(args)), nrx, ncx,
829 
			   COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));
18753 ripley 830 
	else {
831 
	    if(sym)
832 
		symcrossprod(REAL(CAR(args)), nrx, ncx, REAL(ans));
833 
	    else
29387 ripley 834 
		crossprod(REAL(CAR(args)), nrx, ncx,
835 
			  REAL(CADR(args)), nry, ncy, REAL(ans));
18753 ripley 836 
	}
26472 ihaka 837 
 
1820 ihaka 838 
	PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));
26472 ihaka 839 
	if (sym)
18775 ripley 840 
	    PROTECT(ydims = xdims);
841 
	else
842 
	    PROTECT(ydims = getAttrib(CADR(args), R_DimNamesSymbol));
26472 ihaka 843 
 
1820 ihaka 844 
	if (xdims != R_NilValue || ydims != R_NilValue) {
7081 pd 845 
	    SEXP dimnames, dimnamesnames, dnx=R_NilValue, dny=R_NilValue;
26472 ihaka 846 
 
847 
	    /* allocate dimnames and dimnamesnames */
848 
 
6994 pd 849 
	    PROTECT(dimnames = allocVector(VECSXP, 2));
850 
	    PROTECT(dimnamesnames = allocVector(STRSXP, 2));
26472 ihaka 851 
 
6994 pd 852 
	    if (xdims != R_NilValue) {
32350 maechler 853 
		if (ldx == 2) {/* not nrx==1 : .. fixed, ihaka 2003-09-30 */
854 
		    SET_VECTOR_ELT(dimnames, 0, VECTOR_ELT(xdims, 1));
26472 ihaka 855 
		    dnx = getAttrib(xdims, R_NamesSymbol);
856 
		    if(!isNull(dnx))
857 
			SET_STRING_ELT(dimnamesnames, 0, STRING_ELT(dnx, 1));
858 
		}
6994 pd 859 
	    }
26472 ihaka 860 
 
32350 maechler 861 
	    YDIMS_ET_CETERA;
2 r 862 
	}
36679 ripley 863 
 
1820 ihaka 864 
    }
36679 ripley 865 
    else {					/* op == 2: tcrossprod() */
866 
 
867 
	PROTECT(ans = allocMatrix(mode, nrx, nry));
868 
	if (mode == CPLXSXP)
869 
	    if(sym)
870 
		tccrossprod(COMPLEX(CAR(args)), nrx, ncx,
871 
			    COMPLEX(CAR(args)), nry, ncy, COMPLEX(ans));
872 
	    else
873 
		tccrossprod(COMPLEX(CAR(args)), nrx, ncx,
874 
			    COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));
875 
	else {
876 
	    if(sym)
877 
		symtcrossprod(REAL(CAR(args)), nrx, ncx, REAL(ans));
878 
	    else
879 
		tcrossprod(REAL(CAR(args)), nrx, ncx,
880 
			   REAL(CADR(args)), nry, ncy, REAL(ans));
881 
	}
882 
 
883 
	PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));
884 
	if (sym)
885 
	    PROTECT(ydims = xdims);
886 
	else
887 
	    PROTECT(ydims = getAttrib(CADR(args), R_DimNamesSymbol));
888 
 
889 
	if (xdims != R_NilValue || ydims != R_NilValue) {
890 
	    SEXP dimnames, dimnamesnames, dnx=R_NilValue, dny=R_NilValue;
891 
 
892 
	    /* allocate dimnames and dimnamesnames */
893 
 
894 
	    PROTECT(dimnames = allocVector(VECSXP, 2));
895 
	    PROTECT(dimnamesnames = allocVector(STRSXP, 2));
896 
 
897 
	    if (xdims != R_NilValue) {
898 
		if (ldx == 2) {
899 
		    SET_VECTOR_ELT(dimnames, 0, VECTOR_ELT(xdims, 0));
900 
		    dnx = getAttrib(xdims, R_NamesSymbol);
901 
		    if(!isNull(dnx))
902 
			SET_STRING_ELT(dimnamesnames, 0, STRING_ELT(dnx, 0));
903 
		}
904 
	    }
905 
	    if (ydims != R_NilValue) {
906 
		if (ldy == 2) {
907 
		    SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 0));
908 
		    dny = getAttrib(ydims, R_NamesSymbol);
909 
		    if(!isNull(dny))
910 
			SET_STRING_ELT(dimnamesnames, 1, STRING_ELT(dny, 0));
911 
		}
912 
	    }
913 
	    if (VECTOR_ELT(dimnames,0) != R_NilValue ||
914 
		VECTOR_ELT(dimnames,1) != R_NilValue) {
915 
		if (dnx != R_NilValue || dny != R_NilValue)
916 
		    setAttrib(dimnames, R_NamesSymbol, dimnamesnames);
917 
		setAttrib(ans, R_DimNamesSymbol, dimnames);
918 
	    }
919 
 
920 
	    UNPROTECT(2);
921 
	}
922 
    }
1820 ihaka 923 
    UNPROTECT(3);
924 
    return ans;
2 r 925 
}
32350 maechler 926 
#undef YDIMS_ET_CETERA
2 r 927 
 
36990 ripley 928 
SEXP attribute_hidden do_transpose(SEXP call, SEXP op, SEXP args, SEXP rho)
2 r 929 
{
59052 ripley 930 
    SEXP a, r, dims, dimnames, dimnamesnames = R_NilValue,
6994 pd 931 
	ndimnamesnames, rnames, cnames;
59052 ripley 932 
    int ldim, ncol = 0, nrow = 0;
933 
    R_xlen_t len = 0;
2 r 934 
 
1820 ihaka 935 
    checkArity(op, args);
936 
    a = CAR(args);
2 r 937 
 
1839 ihaka 938 
    if (isVector(a)) {
1820 ihaka 939 
	dims = getAttrib(a, R_DimSymbol);
32350 maechler 940 
	ldim = length(dims);
1858 ihaka 941 
	rnames = R_NilValue;
942 
	cnames = R_NilValue;
32350 maechler 943 
	switch(ldim) {
1820 ihaka 944 
	case 0:
59086 ripley 945 
	    len = nrow = LENGTH(a);
1858 ihaka 946 
	    ncol = 1;
947 
	    rnames = getAttrib(a, R_NamesSymbol);
36597 maechler 948 
	    dimnames = rnames;/* for isNull() below*/
1858 ihaka 949 
	    break;
1820 ihaka 950 
	case 1:
59086 ripley 951 
	    len = nrow = LENGTH(a);
1820 ihaka 952 
	    ncol = 1;
32350 maechler 953 
	    dimnames = getAttrib(a, R_DimNamesSymbol);
954 
	    if (dimnames != R_NilValue) {
955 
		rnames = VECTOR_ELT(dimnames, 0);
956 
		dimnamesnames = getAttrib(dimnames, R_NamesSymbol);
957 
	    }
1820 ihaka 958 
	    break;
959 
	case 2:
960 
	    ncol = ncols(a);
961 
	    nrow = nrows(a);
59052 ripley 962 
	    len = XLENGTH(a);
1858 ihaka 963 
	    dimnames = getAttrib(a, R_DimNamesSymbol);
964 
	    if (dimnames != R_NilValue) {
10172 luke 965 
		rnames = VECTOR_ELT(dimnames, 0);
966 
		cnames = VECTOR_ELT(dimnames, 1);
6994 pd 967 
		dimnamesnames = getAttrib(dimnames, R_NamesSymbol);
1858 ihaka 968 
	    }
1820 ihaka 969 
	    break;
970 
	default:
971 
	    goto not_matrix;
2 r 972 
	}
1820 ihaka 973 
    }
1858 ihaka 974 
    else
975 
	goto not_matrix;
1820 ihaka 976 
    PROTECT(r = allocVector(TYPEOF(a), len));
59052 ripley 977 
    R_xlen_t i, j, l_1 = len-1;
1820 ihaka 978 
    switch (TYPEOF(a)) {
979 
    case LGLSXP:
980 
    case INTSXP:
52876 maechler 981 
	// filling in columnwise, "accessing row-wise":
982 
        for (i = 0, j = 0; i < len; i++, j += nrow) {
983 
            if (j > l_1) j -= l_1;
984 
            INTEGER(r)[i] = INTEGER(a)[j];
985 
        }
986 
        break;
1820 ihaka 987 
    case REALSXP:
52876 maechler 988 
        for (i = 0, j = 0; i < len; i++, j += nrow) {
989 
            if (j > l_1) j -= l_1;
990 
            REAL(r)[i] = REAL(a)[j];
991 
        }
992 
        break;
1820 ihaka 993 
    case CPLXSXP:
52876 maechler 994 
        for (i = 0, j = 0; i < len; i++, j += nrow) {
995 
            if (j > l_1) j -= l_1;
996 
            COMPLEX(r)[i] = COMPLEX(a)[j];
997 
        }
998 
        break;
1820 ihaka 999 
    case STRSXP:
52876 maechler 1000 
        for (i = 0, j = 0; i < len; i++, j += nrow) {
1001 
            if (j > l_1) j -= l_1;
1002 
            SET_STRING_ELT(r, i, STRING_ELT(a,j));
1003 
        }
1004 
        break;
1858 ihaka 1005 
    case VECSXP:
52876 maechler 1006 
        for (i = 0, j = 0; i < len; i++, j += nrow) {
1007 
            if (j > l_1) j -= l_1;
1008 
            SET_VECTOR_ELT(r, i, VECTOR_ELT(a,j));
1009 
        }
1010 
        break;
29902 ripley 1011 
    case RAWSXP:
52876 maechler 1012 
        for (i = 0, j = 0; i < len; i++, j += nrow) {
1013 
            if (j > l_1) j -= l_1;
1014 
            RAW(r)[i] = RAW(a)[j];
1015 
        }
1016 
        break;
1858 ihaka 1017 
    default:
52876 maechler 1018 
        UNPROTECT(1);
1019 
        goto not_matrix;
1820 ihaka 1020 
    }
1858 ihaka 1021 
    PROTECT(dims = allocVector(INTSXP, 2));
1820 ihaka 1022 
    INTEGER(dims)[0] = ncol;
1023 
    INTEGER(dims)[1] = nrow;
1024 
    setAttrib(r, R_DimSymbol, dims);
1858 ihaka 1025 
    UNPROTECT(1);
36597 maechler 1026 
    /* R <= 2.2.0: dropped list(NULL,NULL) dimnames :
1027 
     * if(rnames != R_NilValue || cnames != R_NilValue) */
1028 
    if(!isNull(dimnames)) {
1858 ihaka 1029 
	PROTECT(dimnames = allocVector(VECSXP, 2));
10172 luke 1030 
	SET_VECTOR_ELT(dimnames, 0, cnames);
1031 
	SET_VECTOR_ELT(dimnames, 1, rnames);
6994 pd 1032 
	if(!isNull(dimnamesnames)) {
1033 
	    PROTECT(ndimnamesnames = allocVector(VECSXP, 2));
38530 ripley 1034 
	    SET_VECTOR_ELT(ndimnamesnames, 1, STRING_ELT(dimnamesnames, 0));
1035 
	    SET_VECTOR_ELT(ndimnamesnames, 0,
32350 maechler 1036 
			   (ldim == 2) ? STRING_ELT(dimnamesnames, 1):
1037 
			   R_BlankString);
6994 pd 1038 
	    setAttrib(dimnames, R_NamesSymbol, ndimnamesnames);
1039 
	    UNPROTECT(1);
1040 
	}
1858 ihaka 1041 
	setAttrib(r, R_DimNamesSymbol, dimnames);
2 r 1042 
	UNPROTECT(1);
1820 ihaka 1043 
    }
1858 ihaka 1044 
    copyMostAttrib(a, r);
1820 ihaka 1045 
    UNPROTECT(1);
1046 
    return r;
1047 
 not_matrix:
41686 ripley 1048 
    error(_("argument is not a matrix"));
1820 ihaka 1049 
    return call;/* never used; just for -Wall */
2 r 1050 
}
1051 
 
13214 maechler 1052 
/*
1053 
 New version of aperm, using strides for speed.
1054 
 Jonathan Rougier <J.C.Rougier@durham.ac.uk>
1055 
 
1056 
 v1.0 30.01.01
1057 
 
1058 
 M.Maechler : expanded	all ../include/Rdefines.h macros
1059 
 */
1060 
 
1061 
/* this increments iip and sets j using strides */
1062 
 
1063 
#define CLICKJ						\
52649 ripley 1064 
    for (itmp = 0; itmp < n; itmp++)			\
1065 
	if (iip[itmp] == isr[itmp]-1) iip[itmp] = 0;	\
13214 maechler 1066 
	else {						\
1067 
	    iip[itmp]++;				\
1068 
	    break;					\
1069 
	}						\
60305 ripley 1070 
    for (lj = 0, itmp = 0; itmp < n; itmp++)	       	\
1071 
	lj += iip[itmp] * stride[itmp];
17515 ripley 1072 
 
13214 maechler 1073 
/* aperm (a, perm, resize = TRUE) */
36990 ripley 1074 
SEXP attribute_hidden do_aperm(SEXP call, SEXP op, SEXP args, SEXP rho)
2 r 1075 
{
52649 ripley 1076 
    SEXP a, perm, r, dimsa, dimsr, dna;
60305 ripley 1077 
    int i, j, n, itmp;
2 r 1078 
 
1820 ihaka 1079 
    checkArity(op, args);
2 r 1080 
 
1820 ihaka 1081 
    a = CAR(args);
13214 maechler 1082 
    if (!isArray(a))
41686 ripley 1083 
	error(_("invalid first argument, must be an array"));
13214 maechler 1084 
 
1820 ihaka 1085 
    PROTECT(dimsa = getAttrib(a, R_DimSymbol));
13214 maechler 1086 
    n = LENGTH(dimsa);
52649 ripley 1087 
    int *isa = INTEGER(dimsa);
2 r 1088 
 
13214 maechler 1089 
    /* check the permutation */
1090 
 
55420 luke 1091 
    int *pp = (int *) R_alloc((size_t) n, sizeof(int));
52649 ripley 1092 
    perm = CADR(args);
13214 maechler 1093 
    if (length(perm) == 0) {
52649 ripley 1094 
	for (i = 0; i < n; i++) pp[i] = n-1-i;
1095 
    } else {
61696 maechler 1096 
	if (LENGTH(perm) != n)
1097 
	    error(_("'perm' is of wrong length %d (!= %d)"),
1098 
		  LENGTH(perm), n);
1099 
	if (isString(perm)) {
1100 
	    SEXP dna = getAttrib(a, R_DimNamesSymbol);
1101 
	    if (isNull(dna))
1102 
		error(_("'a' does not have named dimnames"));
1103 
	    SEXP dnna = getAttrib(dna, R_NamesSymbol);
1104 
	    if (isNull(dnna))
1105 
		error(_("'a' does not have named dimnames"));
1106 
	    for (i = 0; i < n; i++) {
1107 
		const char *this = translateChar(STRING_ELT(perm, i));
1108 
		for (j = 0; j < n; j++)
1109 
		    if (streql(translateChar(STRING_ELT(dnna, j)),
1110 
			       this)) {pp[i] = j; break;}
1111 
		if (j >= n)
1112 
		    error(_("'perm[%d]' does not match a dimension name"), i+1);
1113 
	    }
1114 
	} else {
1115 
	    PROTECT(perm = coerceVector(perm, INTSXP));
52649 ripley 1116 
	    for (i = 0; i < n; i++) pp[i] = INTEGER(perm)[i] - 1;
1117 
	    UNPROTECT(1);
61696 maechler 1118 
	}
52649 ripley 1119 
    }
2 r 1120 
 
60305 ripley 1121 
    R_xlen_t *iip = (R_xlen_t *) R_alloc((size_t) n, sizeof(R_xlen_t));
41686 ripley 1122 
    for (i = 0; i < n; iip[i++] = 0);
1123 
    for (i = 0; i < n; i++)
52649 ripley 1124 
	if (pp[i] >= 0 && pp[i] < n) iip[pp[i]]++;
1125 
	else error(_("value out of range in 'perm'"));
41686 ripley 1126 
    for (i = 0; i < n; i++)
52649 ripley 1127 
	if (iip[i] == 0) error(_("invalid '%s' argument"), "perm");
2 r 1128 
 
13214 maechler 1129 
    /* create the stride object and permute */
2 r 1130 
 
60305 ripley 1131 
    R_xlen_t *stride = (R_xlen_t *) R_alloc((size_t) n, sizeof(R_xlen_t));
52649 ripley 1132 
    for (iip[0] = 1, i = 1; i<n; i++) iip[i] = iip[i-1] * isa[i-1];
1133 
    for (i = 0; i < n; i++) stride[i] = iip[pp[i]];
13214 maechler 1134 
 
1135 
    /* also need to have the dimensions of r */
1136 
 
52649 ripley 1137 
    PROTECT(dimsr = allocVector(INTSXP, n));
1138 
    int *isr = INTEGER(dimsr);
1139 
    for (i = 0; i < n; i++) isr[i] = isa[pp[i]];
13214 maechler 1140 
 
1141 
    /* and away we go! iip will hold the incrementer */
1142 
 
60305 ripley 1143 
    R_xlen_t len = XLENGTH(a);
1820 ihaka 1144 
    PROTECT(r = allocVector(TYPEOF(a), len));
2 r 1145 
 
52649 ripley 1146 
    for (i = 0; i < n; iip[i++] = 0);
13214 maechler 1147 
 
60305 ripley 1148 
    R_xlen_t li, lj;
1820 ihaka 1149 
    switch (TYPEOF(a)) {
13214 maechler 1150 
 
1820 ihaka 1151 
    case INTSXP:
60305 ripley 1152 
	for (lj = 0, li = 0; li < len; li++) {
1153 
	    INTEGER(r)[li] = INTEGER(a)[lj];
45446 ripley 1154 
	    CLICKJ;
1155 
	}
1156 
	break;
38546 ripley 1157 
 
1820 ihaka 1158 
    case LGLSXP:
60305 ripley 1159 
	for (lj = 0, li = 0; li < len; li++) {
1160 
	    LOGICAL(r)[li] = LOGICAL(a)[lj];
13214 maechler 1161 
	    CLICKJ;
2 r 1162 
	}
1820 ihaka 1163 
	break;
13214 maechler 1164 
 
1820 ihaka 1165 
    case REALSXP:
60305 ripley 1166 
	for (lj = 0, li = 0; li < len; li++) {
1167 
	    REAL(r)[li] = REAL(a)[lj];
13214 maechler 1168 
	    CLICKJ;
1820 ihaka 1169 
	}
1170 
	break;
13214 maechler 1171 
 
1820 ihaka 1172 
    case CPLXSXP:
60305 ripley 1173 
	for (lj = 0, li = 0; li < len; li++) {
1174 
	    COMPLEX(r)[li].r = COMPLEX(a)[lj].r;
1175 
	    COMPLEX(r)[li].i = COMPLEX(a)[lj].i;
13214 maechler 1176 
	    CLICKJ;
1820 ihaka 1177 
	}
1178 
	break;
13214 maechler 1179 
 
1820 ihaka 1180 
    case STRSXP:
60305 ripley 1181 
	for (lj = 0, li = 0; li < len; li++) {
1182 
	    SET_STRING_ELT(r, li, STRING_ELT(a, lj));
13214 maechler 1183 
	    CLICKJ;
1820 ihaka 1184 
	}
1185 
	break;
13214 maechler 1186 
 
1839 ihaka 1187 
    case VECSXP:
60305 ripley 1188 
	for (lj = 0, li = 0; li < len; li++) {
1189 
	    SET_VECTOR_ELT(r, li, VECTOR_ELT(a, lj));
13214 maechler 1190 
	    CLICKJ;
1839 ihaka 1191 
	}
12976 pd 1192 
	break;
13214 maechler 1193 
 
29902 ripley 1194 
    case RAWSXP:
60305 ripley 1195 
	for (lj = 0, li = 0; li < len; li++) {
1196 
	    RAW(r)[li] = RAW(a)[lj];
29902 ripley 1197 
	    CLICKJ;
1198 
	}
1199 
	break;
1200 
 
1820 ihaka 1201 
    default:
31908 ripley 1202 
	UNIMPLEMENTED_TYPE("aperm", a);
1820 ihaka 1203 
    }
17515 ripley 1204 
 
13214 maechler 1205 
    /* handle the resize */
52649 ripley 1206 
    int resize = asLogical(CADDR(args));
1207 
    if (resize == NA_LOGICAL) error(_("'resize' must be TRUE or FALSE"));
1208 
    setAttrib(r, R_DimSymbol, resize ? dimsr : dimsa);
13214 maechler 1209 
 
52649 ripley 1210 
    /* and handle the dimnames, if any */
1211 
    if (resize) {
1212 
	PROTECT(dna = getAttrib(a, R_DimNamesSymbol));
1213 
	if (dna != R_NilValue) {
1214 
	    SEXP dnna, dnr, dnnr;
13214 maechler 1215 
 
52649 ripley 1216 
	    PROTECT(dnr  = allocVector(VECSXP, n));
1217 
	    PROTECT(dnna = getAttrib(dna, R_NamesSymbol));
1218 
	    if (dnna != R_NilValue) {
1219 
		PROTECT(dnnr = allocVector(STRSXP, n));
1220 
		for (i = 0; i < n; i++) {
1221 
		    SET_VECTOR_ELT(dnr, i, VECTOR_ELT(dna, pp[i]));
1222 
		    SET_STRING_ELT(dnnr, i, STRING_ELT(dnna, pp[i]));
1223 
		}
1224 
		setAttrib(dnr, R_NamesSymbol, dnnr);
1225 
		UNPROTECT(1);
1226 
	    } else {
1227 
		for (i = 0; i < n; i++)
1228 
		    SET_VECTOR_ELT(dnr, i, VECTOR_ELT(dna, pp[i]));
1229 
	    }
1230 
	    setAttrib(r, R_DimNamesSymbol, dnr);
1231 
	    UNPROTECT(2);
13214 maechler 1232 
	}
52649 ripley 1233 
	UNPROTECT(1);
13214 maechler 1234 
    }
17515 ripley 1235 
 
52649 ripley 1236 
    UNPROTECT(3); /* dimsa, r, dimsr */
1820 ihaka 1237 
    return r;
2 r 1238 
}
17515 ripley 1239 
 
1240 
/* colSums(x, n, p, na.rm) and friends */
36990 ripley 1241 
SEXP attribute_hidden do_colsum(SEXP call, SEXP op, SEXP args, SEXP rho)
17515 ripley 1242 
{
1243 
    SEXP x, ans = R_NilValue;
60282 ripley 1244 
    int type;
17521 ripley 1245 
    Rboolean NaRm, keepNA;
17515 ripley 1246 
 
1247 
    checkArity(op, args);
1248 
    x = CAR(args); args = CDR(args);
60282 ripley 1249 
    int n = asInteger(CAR(args)); args = CDR(args);
1250 
    int p = asInteger(CAR(args)); args = CDR(args);
17515 ripley 1251 
    NaRm = asLogical(CAR(args));
33840 ripley 1252 
    if (n == NA_INTEGER || n < 0)
44512 ripley 1253 
	error(_("invalid '%s' argument"), "n");
33840 ripley 1254 
    if (p == NA_INTEGER || p < 0)
44512 ripley 1255 
	error(_("invalid '%s' argument"), "p");
1256 
    if (NaRm == NA_LOGICAL) error(_("invalid '%s' argument"), "na.rm");
17521 ripley 1257 
    keepNA = !NaRm;
17515 ripley 1258 
 
60282 ripley 1259 
    int OP = PRIMVAL(op);
17515 ripley 1260 
    switch (type = TYPEOF(x)) {
1261 
    case LGLSXP: break;
1262 
    case INTSXP: break;
1263 
    case REALSXP: break;
1264 
    default:
41686 ripley 1265 
	error(_("'x' must be numeric"));
17515 ripley 1266 
    }
1267 
 
17521 ripley 1268 
    if (OP == 0 || OP == 1) { /* columns */
17515 ripley 1269 
	PROTECT(ans = allocVector(REALSXP, p));
62637 ripley 1270 
#ifdef _OPENMP
60282 ripley 1271 
	int nthreads;
55099 ripley 1272 
	/* This gives a spurious -Wunused-but-set-variable error */
53854 luke 1273 
	if (R_num_math_threads > 0)
1274 
	    nthreads = R_num_math_threads;
1275 
	else
1276 
	    nthreads = 1; /* for now */
1277 
#pragma omp parallel for num_threads(nthreads) default(none) \
59537 luke 1278 
    firstprivate(x, ans, n, p, type, NaRm, keepNA, R_NaReal, R_NaInt, OP)
53854 luke 1279 
#endif
59537 luke 1280 
	for (int j = 0; j < p; j++) {
1281 
	    int cnt = n, i;
60268 ripley 1282 
	    LDOUBLE sum = 0.0;
17515 ripley 1283 
	    switch (type) {
1284 
	    case REALSXP:
60282 ripley 1285 
	    {
1286 
		double *rx = REAL(x) + (R_xlen_t)n*j;
17521 ripley 1287 
		if (keepNA)
1288 
		    for (sum = 0., i = 0; i < n; i++) sum += *rx++;
1289 
		else {
1290 
		    for (cnt = 0, sum = 0., i = 0; i < n; i++, rx++)
1291 
			if (!ISNAN(*rx)) {cnt++; sum += *rx;}
60282 ripley 1292 
			else if (keepNA) {sum = NA_REAL; break;} // unused
17521 ripley 1293 
		}
17515 ripley 1294 
		break;
60282 ripley 1295 
	    }
17515 ripley 1296 
	    case INTSXP:
60282 ripley 1297 
	    {
1298 
		int *ix = INTEGER(x) + (R_xlen_t)n*j;
17521 ripley 1299 
		for (cnt = 0, sum = 0., i = 0; i < n; i++, ix++)
1300 
		    if (*ix != NA_INTEGER) {cnt++; sum += *ix;}
1301 
		    else if (keepNA) {sum = NA_REAL; break;}
17515 ripley 1302 
		break;
60282 ripley 1303 
	    }
17515 ripley 1304 
	    case LGLSXP:
60282 ripley 1305 
	    {
1306 
		int *ix = LOGICAL(x) + (R_xlen_t)n*j;
17521 ripley 1307 
		for (cnt = 0, sum = 0., i = 0; i < n; i++, ix++)
1308 
		    if (*ix != NA_LOGICAL) {cnt++; sum += *ix;}
1309 
		    else if (keepNA) {sum = NA_REAL; break;}
17515 ripley 1310 
		break;
1311 
	    }
60282 ripley 1312 
	    }
59537 luke 1313 
	    if (OP == 1) sum /= cnt; /* gives NaN for cnt = 0 */
59095 ripley 1314 
	    REAL(ans)[j] = (double) sum;
17515 ripley 1315 
	}
1316 
    }
59537 luke 1317 
    else { /* rows */
17515 ripley 1318 
	PROTECT(ans = allocVector(REALSXP, n));
17521 ripley 1319 
 
59537 luke 1320 
	/* allocate scratch storage to allow accumulating by columns
1321 
	   to improve cache hits */
1322 
	int *Cnt = NULL;
60268 ripley 1323 
	LDOUBLE *rans;
59537 luke 1324 
	if(n <= 10000) {
60268 ripley 1325 
	    R_CheckStack2(n * sizeof(LDOUBLE));
1326 
	    rans = (LDOUBLE *) alloca(n * sizeof(LDOUBLE));
60517 ripley 1327 
	    Memzero(rans, n);
60268 ripley 1328 
	} else rans = Calloc(n, LDOUBLE);
59537 luke 1329 
	if (!keepNA && OP == 3) Cnt = Calloc(n, int);
1330 
 
60282 ripley 1331 
	for (int j = 0; j < p; j++) {
60268 ripley 1332 
	    LDOUBLE *ra = rans;
59537 luke 1333 
	    switch (type) {
1334 
	    case REALSXP:
60282 ripley 1335 
	    {
1336 
		double *rx = REAL(x) + (R_xlen_t)n * j;
17521 ripley 1337 
		if (keepNA)
59537 luke 1338 
		    for (int i = 0; i < n; i++) *ra++ += *rx++;
17521 ripley 1339 
		else
59537 luke 1340 
		    for (int i = 0; i < n; i++, ra++, rx++)
17521 ripley 1341 
			if (!ISNAN(*rx)) {
1342 
			    *ra += *rx;
59537 luke 1343 
			    if (OP == 3) Cnt[i]++;
17521 ripley 1344 
			}
59537 luke 1345 
		break;
60282 ripley 1346 
	    }
17515 ripley 1347 
	    case INTSXP:
60282 ripley 1348 
	    {
1349 
		int *ix = INTEGER(x) + (R_xlen_t)n * j;
59537 luke 1350 
		for (int i = 0; i < n; i++, ra++, ix++)
1351 
		    if (keepNA) {
1352 
			if (*ix != NA_INTEGER) *ra += *ix;
1353 
			else *ra = NA_REAL;
1354 
		    }
1355 
		    else if (*ix != NA_INTEGER) {
1356 
			*ra += *ix;
1357 
			if (OP == 3) Cnt[i]++;
1358 
		    }
17515 ripley 1359 
		break;
60282 ripley 1360 
	    }
17515 ripley 1361 
	    case LGLSXP:
60282 ripley 1362 
	    {
1363 
		int *ix = LOGICAL(x) + (R_xlen_t)n * j;
59537 luke 1364 
		for (int i = 0; i < n; i++, ra++, ix++)
1365 
		    if (keepNA) {
1366 
			if (*ix != NA_LOGICAL) *ra += *ix;
1367 
			else *ra = NA_REAL;
1368 
		    }
1369 
		    else if (*ix != NA_LOGICAL) {
1370 
			*ra += *ix;
1371 
			if (OP == 3) Cnt[i]++;
1372 
		    }
17515 ripley 1373 
		break;
1374 
	    }
60282 ripley 1375 
	    }
17515 ripley 1376 
	}
59537 luke 1377 
	if (OP == 3) {
1378 
	    if (keepNA)
1379 
		for (int i = 0; i < n; i++) rans[i] /= p;
61696 maechler 1380 
	    else
59537 luke 1381 
		for (int i = 0; i < n; i++) rans[i] /= Cnt[i];
1382 
	}
1383 
	for (int i = 0; i < n; i++) REAL(ans)[i] = (double) rans[i];
1384 
 
1385 
	if (!keepNA && OP == 3) Free(Cnt);
1386 
	if(n > 10000) Free(rans);
17515 ripley 1387 
    }
1388 
 
1389 
    UNPROTECT(1);
1390 
    return ans;
1391 
}
59621 ripley 1392 
 
1393 
/*
1394 
{
1395 
    data <- as.vector(data)
1396 
    dim <- as.integer(dim)
1397 
    vl <- prod(dim)
1398 
    if (length(data) != vl) {
61696 maechler 1399 
        if (vl > .Machine$integer.max)
59621 ripley 1400 
            stop("'dim' specifies too large an array")
1401 
        data <- rep(data, length.out = vl)
1402 
    }
61696 maechler 1403 
    if (length(dim))
59621 ripley 1404 
        dim(data) <- dim
61696 maechler 1405 
    if (is.list(dimnames) && length(dimnames))
59621 ripley 1406 
        dimnames(data) <- dimnames
1407 
    data
1408 
}
1409 
*/
1410 
 
1411 
/* array(data, dim, dimnames) */
1412 
SEXP attribute_hidden do_array(SEXP call, SEXP op, SEXP args, SEXP rho)
1413 
{
1414 
    SEXP vals, ans, dims, dimnames;
1415 
    R_xlen_t lendat, i, nans;
1416 
 
1417 
    checkArity(op, args);
1418 
    vals = CAR(args);
61566 ripley 1419 
    /* at least NULL can get here */
59621 ripley 1420 
    switch(TYPEOF(vals)) {
1421 
	case LGLSXP:
1422 
	case INTSXP:
1423 
	case REALSXP:
1424 
	case CPLXSXP:
1425 
	case STRSXP:
1426 
	case RAWSXP:
1427 
	case EXPRSXP:
1428 
	case VECSXP:
1429 
	    break;
1430 
	default:
61566 ripley 1431 
	    error(_("'data' must be of a vector type, was '%s'"),
1432 
		type2char(TYPEOF(vals)));
59621 ripley 1433 
    }
1434 
    lendat = XLENGTH(vals);
1435 
    dims = CADR(args);
1436 
    dimnames = CADDR(args);
1437 
    PROTECT(dims = coerceVector(dims, INTSXP));
1438 
    int nd = LENGTH(dims);
59622 ripley 1439 
    if (nd == 0) error(_("'dims' cannot be of length 0"));
1440 
    double d = 1.0;
1441 
    for (int j = 0; j < nd; j++) d *= INTEGER(dims)[j];
59621 ripley 1442 
#ifndef LONG_VECTOR_SUPPORT
59622 ripley 1443 
    if (d > INT_MAX) error(_("too many elements specified"));
59621 ripley 1444 
#endif
59622 ripley 1445 
    nans = (R_xlen_t) d;
59621 ripley 1446 
 
1447 
    PROTECT(ans = allocVector(TYPEOF(vals), nans));
1448 
    switch(TYPEOF(vals)) {
1449 
    case LGLSXP:
1450 
	if (nans && lendat)
1451 
	    for (i = 0; i < nans; i++)
1452 
		LOGICAL(ans)[i] = LOGICAL(vals)[i % lendat];
1453 
	else
1454 
	    for (i = 0; i < nans; i++) LOGICAL(ans)[i] = NA_LOGICAL;
1455 
	break;
1456 
    case INTSXP:
1457 
	if (nans && lendat)
1458 
	    for (i = 0; i < nans; i++)
1459 
		INTEGER(ans)[i] = INTEGER(vals)[i % lendat];
1460 
	else
1461 
	    for (i = 0; i < nans; i++) INTEGER(ans)[i] = NA_INTEGER;
1462 
	break;
1463 
    case REALSXP:
1464 
	if (nans && lendat)
1465 
	    for (i = 0; i < nans; i++) REAL(ans)[i] = REAL(vals)[i % lendat];
1466 
	else
1467 
	    for (i = 0; i < nans; i++) REAL(ans)[i] = NA_REAL;
1468 
	break;
1469 
    case CPLXSXP:
1470 
	if (nans && lendat)
1471 
	    for (i = 0; i < nans; i++)
1472 
		COMPLEX(ans)[i] = COMPLEX(vals)[i % lendat];
1473 
	else {
1474 
	    Rcomplex na_cmplx;
1475 
	    na_cmplx.r = NA_REAL;
1476 
	    na_cmplx.i = 0;
1477 
	    for (i = 0; i < nans; i++) COMPLEX(ans)[i] = na_cmplx;
1478 
	}
1479 
	break;
1480 
    case RAWSXP:
1481 
	if (nans && lendat)
1482 
	    for (i = 0; i < nans; i++) RAW(ans)[i] = RAW(vals)[i % lendat];
1483 
	else
1484 
	    for (i = 0; i < nans; i++) RAW(ans)[i] = 0;
1485 
	break;
1486 
    /* Rest are already initialized */
1487 
    case STRSXP:
1488 
	if (nans && lendat)
1489 
	    for (i = 0; i < nans; i++)
1490 
		SET_STRING_ELT(ans, i, STRING_ELT(vals, i % lendat));
1491 
	break;
1492 
    case VECSXP:
1493 
    case EXPRSXP:
1494 
	if (nans && lendat)
1495 
	    for (i = 0; i < nans; i++)
60432 ripley 1496 
		SET_VECTOR_ELT(ans, i, VECTOR_ELT(vals, i % lendat));
59621 ripley 1497 
	break;
1498 
    default:
1499 
	// excluded above
1500 
	break;
1501 
    }
1502 
 
59622 ripley 1503 
    ans = dimgets(ans, dims);
1504 
    if (TYPEOF(dimnames) == VECSXP && LENGTH(dimnames)) {
1505 
	PROTECT(ans);
1506 
	ans = dimnamesgets(ans, dimnames);
1507 
	UNPROTECT(1);
1508 
    }
61696 maechler 1509 
 
59621 ripley 1510 
    UNPROTECT(2);
1511 
    return ans;
1512 
}
60287 ripley 1513 
 
1514 
SEXP attribute_hidden do_diag(SEXP call, SEXP op, SEXP args, SEXP rho)
1515 
{
1516 
    SEXP ans, x, snr, snc;
1517 
    int nr = 1, nc = 1, nprotect = 1;
1518 
 
1519 
    checkArity(op, args);
1520 
    x = CAR(args);
1521 
    snr = CADR(args);
1522 
    snc = CADDR(args);
1523 
    nr = asInteger(snr);
1524 
    if (nr == NA_INTEGER)
1525 
	error(_("invalid 'nrow' value (too large or NA)"));
1526 
    if (nr < 0)
1527 
	error(_("invalid 'nrow' value (< 0)"));
1528 
    nc = asInteger(snc);
1529 
    if (nc == NA_INTEGER)
1530 
	error(_("invalid 'ncol' value (too large or NA)"));
1531 
    if (nc < 0)
1532 
	error(_("invalid 'ncol' value (< 0)"));
61696 maechler 1533 
    int mn = (nr < nc) ? nr : nc;
60287 ripley 1534 
    if (mn > 0 && LENGTH(x) == 0)
1535 
	error(_("'x' must have positive length"));
1536 
 
1537 
 #ifndef LONG_VECTOR_SUPPORT
1538 
   if ((double)nr * (double)nc > INT_MAX)
1539 
	error(_("too many elements specified"));
1540 
#endif
1541 
 
1542 
   if (TYPEOF(x) == CPLXSXP) {
1543 
       PROTECT(ans = allocMatrix(CPLXSXP, nr, nc));
1544 
       int nx = LENGTH(x);
1545 
       R_xlen_t NR = nr;
1546 
       Rcomplex *rx = COMPLEX(x), *ra = COMPLEX(ans), zero;
1547 
       zero.r = zero.i = 0.0;
1548 
       for (R_xlen_t i = 0; i < NR*nc; i++) ra[i] = zero;
1549 
       for (int j = 0; j < mn; j++) ra[j * (NR+1)] = rx[j % nx];
1550 
  } else {
1551 
       if(TYPEOF(x) != REALSXP) {
1552 
	   PROTECT(x = coerceVector(x, REALSXP));
1553 
	   nprotect++;
1554 
       }
1555 
       PROTECT(ans = allocMatrix(REALSXP, nr, nc));
1556 
       int nx = LENGTH(x);
1557 
       R_xlen_t NR = nr;
1558 
       double *rx = REAL(x), *ra = REAL(ans);
1559 
       for (R_xlen_t i = 0; i < NR*nc; i++) ra[i] = 0.0;
1560 
       for (int j = 0; j < mn; j++) ra[j * (NR+1)] = rx[j % nx];
1561 
   }
1562 
   UNPROTECT(nprotect);
1563 
   return ans;
1564 
}
60321 ripley 1565 
 
1566 
 
60335 ripley 1567 
/* backsolve(r, b, k, upper.tri, transpose) */
60321 ripley 1568 
SEXP attribute_hidden do_backsolve(SEXP call, SEXP op, SEXP args, SEXP rho)
1569 
{
1570 
    int nprot = 1;
1571 
    checkArity(op, args);
1572 
 
1573 
    SEXP r = CAR(args); args = CDR(args);
60335 ripley 1574 
    SEXP b = CAR(args); args = CDR(args);
1575 
    int nrr = nrows(r), nrb = nrows(b), ncb = ncols(b);
60321 ripley 1576 
    int k = asInteger(CAR(args)); args = CDR(args);
60335 ripley 1577 
    /* k is the number of rows to be used: there must be at least that
1578 
       many rows and cols in the rhs and at least that many rows on
1579 
       the rhs.
1580 
    */
61696 maechler 1581 
    if (k == NA_INTEGER || k <= 0 || k > nrr || k > ncols(r) || k > nrb)
60321 ripley 1582 
	error(_("invalid '%s' argument"), "k");
1583 
    int upper = asLogical(CAR(args)); args = CDR(args);
1584 
    if (upper == NA_INTEGER) error(_("invalid '%s' argument"), "upper.tri");
1585 
    int trans = asLogical(CAR(args));
1586 
    if (trans == NA_INTEGER) error(_("invalid '%s' argument"), "transpose");
1587 
    if (TYPEOF(r) != REALSXP) {PROTECT(r = coerceVector(r, REALSXP)); nprot++;}
60335 ripley 1588 
    if (TYPEOF(b) != REALSXP) {PROTECT(b = coerceVector(b, REALSXP)); nprot++;}
1589 
    double *rr = REAL(r);
1590 
 
1591 
    /* check for zeros on diagonal of r: only k row/cols are used. */
1592 
    size_t incr = nrr + 1;
1593 
    for(int i = 0; i < k; i++) { /* check for zeros on diagonal */
1594 
	if (rr[i * incr] == 0.0)
1595 
	    error(_("singular matrix in 'backsolve'. First zero in diagonal [%d]"),
1596 
		  i + 1);
1597 
    }
1598 
 
62618 ripley 1599 
    SEXP ans = PROTECT(allocMatrix(REALSXP, k, ncb));
60335 ripley 1600 
    if (k > 0 && ncb > 0) {
1601 
       /* copy (part) cols of b to ans */
1602 
	for(R_xlen_t j = 0; j < ncb; j++)
1603 
	    memcpy(REAL(ans) + j*k, REAL(b) + j*nrb, (size_t)k *sizeof(double));
1604 
	double one = 1.0;
61696 maechler 1605 
	F77_CALL(dtrsm)("L", upper ? "U" : "L", trans ? "T" : "N", "N",
60335 ripley 1606 
			&k, &ncb, &one, rr, &nrr, REAL(ans), &k);
1607 
    }
60321 ripley 1608 
    UNPROTECT(nprot);
1609 
    return ans;
1610 
}
60327 ripley 1611 
 
1612 
/* max.col(m, ties.method) */
1613 
SEXP attribute_hidden do_maxcol(SEXP call, SEXP op, SEXP args, SEXP rho)
1614 
{
1615 
    checkArity(op, args);
1616 
    SEXP m = CAR(args);
1617 
    int method = asInteger(CADR(args));
1618 
    int nr = nrows(m), nc = ncols(m), nprot = 1;
62618 ripley 1619 
    if (TYPEOF(m) != REALSXP) {PROTECT(m = coerceVector(m, REALSXP)); nprot++;}
1620 
    SEXP ans = PROTECT(allocVector(INTSXP, nr));
60327 ripley 1621 
    R_max_col(REAL(m), &nr, &nc, INTEGER(ans), &method);
1622 
    UNPROTECT(nprot);
1623 
    return ans;
1624 
}