The R Project SVN R

Rev

Details | 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
26467 maechler 4 
 *  Copyright (C) 1998-2001   The R Development Core Team
29387 ripley 5 
 *  Copyright (C) 2002--2004  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 
 *
26467 maechler 17 
 *  A copy of the GNU General Public License is available via WWW at
18 
 *  http://www.gnu.org/copyleft/gpl.html.  You can also obtain it by
19 
 *  writing to the Free Software Foundation, Inc., 59 Temple Place,
20 
 *  Suite 330, Boston, MA  02111-1307  USA.
2 r 21 
 */
22 
 
5187 hornik 23 
#ifdef HAVE_CONFIG_H
7701 hornik 24 
#include <config.h>
5187 hornik 25 
#endif
26 
 
11499 ripley 27 
#include <Defn.h>
28 
#include <Rmath.h>
13942 bates 29 
#include <R_ext/RS.h>
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.''
41 
*/
1839 ihaka 42 
SEXP GetRowNames(SEXP dimnames)
43 
{
44 
    if (TYPEOF(dimnames) == VECSXP)
10172 luke 45 
	return VECTOR_ELT(dimnames, 0);
1839 ihaka 46 
    else
1858 ihaka 47 
	return R_NilValue;
2213 maechler 48 
}
1839 ihaka 49 
 
50 
SEXP GetColNames(SEXP dimnames)
51 
{
52 
    if (TYPEOF(dimnames) == VECSXP)
10172 luke 53 
	return VECTOR_ELT(dimnames, 1);
1839 ihaka 54 
    else
2213 maechler 55 
	return R_NilValue;
1839 ihaka 56 
}
57 
 
2 r 58 
SEXP do_matrix(SEXP call, SEXP op, SEXP args, SEXP rho)
59 
{
1820 ihaka 60 
    SEXP vals, snr, snc;
61 
    int nr, nc, byrow, lendat;
2 r 62 
 
1820 ihaka 63 
    checkArity(op, args);
64 
    vals = CAR(args);
65 
    snr = CADR(args);
66 
    snc = CADDR(args);
67 
    byrow = asInteger(CADR(CDDR(args)));
2 r 68 
 
28014 ripley 69 
    /* R wrapper does as.vector
1820 ihaka 70 
    if (isVector(vals) || isList(vals)) {
28014 ripley 71 
	if (length(vals) < 0)  (sic! cannot happen)
72 
	   errorcall(call, "argument has length zero");
1881 ihaka 73 
    }
28014 ripley 74 
    else errorcall(call, "invalid matrix element type"); */
2 r 75 
 
1820 ihaka 76 
    if (!isNumeric(snr) || !isNumeric(snc))
5731 ripley 77 
	error("non-numeric matrix extent");
2 r 78 
 
1820 ihaka 79 
    lendat = length(vals);
80 
    nr = asInteger(snr);
28349 ripley 81 
    if (nr == NA_INTEGER) /* This is < 0 */
82 
	error("matrix: invalid nrow value (too large or NA)");
83 
    if (nr < 0)
84 
	error("matrix: invalid nrow value (< 0)");
1820 ihaka 85 
    nc = asInteger(snc);
28349 ripley 86 
    if (nc < 0)
87 
	error("matrix: invalid ncol value (< 0)");
88 
    if (nc == NA_INTEGER)
89 
	error("matrix: invalid ncol value (too large or NA)");
90 
    if (nc < 0)
91 
	error("matrix: invalid ncol value (< 0)");
2 r 92 
 
28014 ripley 93 
    if(lendat > 0 ) {
94 
	if (lendat > 1 && (nr * nc) % lendat != 0) {
95 
	    if (((lendat > nr) && (lendat / nr) * nr != lendat) ||
96 
		((lendat < nr) && (nr / lendat) * lendat != nr))
97 
		warning("data length [%d] is not a sub-multiple or multiple of the number of rows [%d] in matrix", lendat, nr);
98 
	    else if (((lendat > nc) && (lendat / nc) * nc != lendat) ||
99 
		     ((lendat < nc) && (nc / lendat) * lendat != nc))
100 
		warning("data length [%d] is not a sub-multiple or multiple of the number of columns [%d] in matrix", lendat, nc);
101 
	}
1839 ihaka 102 
	else if ((lendat > 1) && (nr * nc == 0)){
28014 ripley 103 
	    warning("data length exceeds size of matrix");
2 r 104 
	}
28014 ripley 105 
    }
2 r 106 
 
28349 ripley 107 
    if ((double)nr * (double)nc > INT_MAX)
108 
	error("matrix: too many elements specified");
109 
 
1820 ihaka 110 
    PROTECT(snr = allocMatrix(TYPEOF(vals), nr, nc));
28014 ripley 111 
    if(lendat) {
112 
	if (isVector(vals))
113 
	    copyMatrix(snr, vals, byrow);
114 
	else
115 
	    copyListMatrix(snr, vals, byrow);
116 
    } else if (isVector(vals)) { /* fill with NAs */
117 
	int i, j;
118 
	switch(TYPEOF(vals)) {
119 
	case STRSXP:
120 
	    for (i = 0; i < nr; i++)
121 
		for (j = 0; j < nc; j++)
122 
		    SET_STRING_ELT(snr, i + j * nr, NA_STRING);
123 
	    break;
124 
	case LGLSXP:
125 
	    for (i = 0; i < nr; i++)
126 
		for (j = 0; j < nc; j++)
127 
		    LOGICAL(snr)[i + j * nr] = NA_LOGICAL;
128 
	    break;
129 
	case INTSXP:
130 
	    for (i = 0; i < nr; i++)
131 
		for (j = 0; j < nc; j++)
132 
		    INTEGER(snr)[i + j * nr] = NA_INTEGER;
133 
	    break;
134 
	case REALSXP:
135 
	    for (i = 0; i < nr; i++)
136 
		for (j = 0; j < nc; j++)
137 
		    REAL(snr)[i + j * nr] = NA_REAL;
138 
	    break;
139 
	case CPLXSXP:
140 
	    {
141 
		Rcomplex na_cmplx;
142 
		na_cmplx.r = NA_REAL;
143 
		na_cmplx.i = 0;
144 
		for (i = 0; i < nr; i++)
145 
		    for (j = 0; j < nc; j++)
146 
			COMPLEX(snr)[i + j * nr] = na_cmplx;
147 
	    }
148 
	    break;
30691 murdoch 149 
	case RAWSXP:
150 
	    for (i = 0; i < nr; i++)
151 
		for (j = 0; j < nc; j++)
152 
		    RAW(snr)[i + j * nr] = 0;
153 
	    break;
28014 ripley 154 
	}
155 
    }
1820 ihaka 156 
    UNPROTECT(1);
157 
    return snr;
2 r 158 
}
159 
 
160 
 
161 
SEXP allocMatrix(SEXPTYPE mode, int nrow, int ncol)
162 
{
1820 ihaka 163 
    SEXP s, t;
164 
    int n;
2 r 165 
 
1820 ihaka 166 
    if (nrow < 0 || ncol < 0)
5731 ripley 167 
	error("negative extents to matrix");
28349 ripley 168 
    if ((double)nrow * (double)ncol > INT_MAX)
169 
	error("allocMatrix: too many elements specified");
1820 ihaka 170 
    n = nrow * ncol;
171 
    PROTECT(s = allocVector(mode, n));
172 
    PROTECT(t = allocVector(INTSXP, 2));
173 
    INTEGER(t)[0] = nrow;
174 
    INTEGER(t)[1] = ncol;
175 
    setAttrib(s, R_DimSymbol, t);
176 
    UNPROTECT(2);
177 
    return s;
2 r 178 
}
179 
 
1130 maechler 180 
 
2 r 181 
SEXP allocArray(SEXPTYPE mode, SEXP dims)
182 
{
1820 ihaka 183 
    SEXP array;
184 
    int i, n;
28349 ripley 185 
    double dn;
2 r 186 
 
28349 ripley 187 
    dn = n = 1;
188 
    for (i = 0; i < LENGTH(dims); i++) {
189 
	dn *= INTEGER(dims)[i];
190 
	if(dn > INT_MAX)
191 
	    error("allocArray: too many elements specified by dims");
192 
	n *= INTEGER(dims)[i];
193 
    }
2 r 194 
 
1820 ihaka 195 
    PROTECT(dims = duplicate(dims));
196 
    PROTECT(array = allocVector(mode, n));
197 
    setAttrib(array, R_DimSymbol, dims);
198 
    UNPROTECT(2);
199 
    return array;
2 r 200 
}
201 
 
1820 ihaka 202 
/* DropDims strips away redundant dimensioning information. */
203 
/* If there is an appropriate dimnames attribute the correct */
204 
/* element is extracted and attached to the vector as a names */
205 
/* attribute.  Note that this function mutates x. */
206 
/* Duplication should occur before this is called. */
2 r 207 
 
208 
SEXP DropDims(SEXP x)
209 
{
2506 maechler 210 
    SEXP q, dims, dimnames, newnames = R_NilValue;
1839 ihaka 211 
    int i, n, ndims;
2 r 212 
 
1820 ihaka 213 
    PROTECT(x);
214 
    dims = getAttrib(x, R_DimSymbol);
215 
    dimnames = getAttrib(x, R_DimNamesSymbol);
2 r 216 
 
1820 ihaka 217 
    /* Check that dropping will actually do something. */
218 
    /* (1) Check that there is a "dim" attribute. */
2 r 219 
 
1839 ihaka 220 
    if (dims == R_NilValue) {
1820 ihaka 221 
	UNPROTECT(1);
222 
	return x;
223 
    }
1839 ihaka 224 
    ndims = LENGTH(dims);
2 r 225 
 
1839 ihaka 226 
    /* (2) Check whether there are redundant extents */
1820 ihaka 227 
    n = 0;
1839 ihaka 228 
    for (i = 0; i < ndims; i++)
229 
	if (INTEGER(dims)[i] != 1) n++;
230 
    if (n == ndims) {
1820 ihaka 231 
	UNPROTECT(1);
232 
	return x;
233 
    }
2 r 234 
 
1839 ihaka 235 
    if (n <= 1) {
236 
	/* We have reduced to a vector result. */
1820 ihaka 237 
	if (dimnames != R_NilValue) {
238 
	    n = length(dims);
1839 ihaka 239 
	    if (TYPEOF(dimnames) == VECSXP) {
240 
		for (i = 0; i < n; i++) {
241 
		    if (INTEGER(dims)[i] != 1) {
10172 luke 242 
			newnames = VECTOR_ELT(dimnames, i);
1839 ihaka 243 
			break;
244 
		    }
2 r 245 
		}
1820 ihaka 246 
	    }
1839 ihaka 247 
	    else {
248 
		q = dimnames;
249 
		for (i = 0; i < n; i++) {
250 
		    if (INTEGER(dims)[i] != 1) {
251 
			newnames = CAR(q);
252 
			break;
253 
		    }
254 
		    q = CDR(q);
1820 ihaka 255 
		}
256 
	    }
2 r 257 
	}
1820 ihaka 258 
	PROTECT(newnames);
259 
	setAttrib(x, R_DimNamesSymbol, R_NilValue);
260 
	setAttrib(x, R_DimSymbol, R_NilValue);
261 
	setAttrib(x, R_NamesSymbol, newnames);
262 
	UNPROTECT(1);
263 
    }
1839 ihaka 264 
    else {
265 
	/* We have a lower dimensional array. */
6994 pd 266 
	SEXP newdims, dnn, newnamesnames = R_NilValue;
267 
	dnn = getAttrib(dimnames, R_NamesSymbol);
1820 ihaka 268 
	PROTECT(newdims = allocVector(INTSXP, n));
1858 ihaka 269 
	for (i = 0, n = 0; i < ndims; i++)
1839 ihaka 270 
	    if (INTEGER(dims)[i] != 1)
1820 ihaka 271 
		INTEGER(newdims)[n++] = INTEGER(dims)[i];
6994 pd 272 
	if (!isNull(dimnames)) {
1858 ihaka 273 
	    int havenames = 0;
274 
	    for (i = 0; i < ndims; i++)
10172 luke 275 
		if (INTEGER(dims)[i] != 1 &&
276 
		    VECTOR_ELT(dimnames, i) != R_NilValue)
1858 ihaka 277 
		    havenames = 1;
278 
	    if (havenames) {
2715 pd 279 
		PROTECT(newnames = allocVector(VECSXP, n));
6994 pd 280 
		PROTECT(newnamesnames = allocVector(STRSXP, n));
281 
		for (i = 0, n = 0; i < ndims; i++) {
282 
		    if (INTEGER(dims)[i] != 1) {
283 
			if(!isNull(dnn))
10172 luke 284 
			    SET_STRING_ELT(newnamesnames, n,
285 
					   STRING_ELT(dnn, i));
286 
			SET_VECTOR_ELT(newnames, n++, VECTOR_ELT(dimnames, i));
6994 pd 287 
		    }
1858 ihaka 288 
		}
1820 ihaka 289 
	    }
1858 ihaka 290 
	    else dimnames = R_NilValue;
1839 ihaka 291 
	}
1858 ihaka 292 
	PROTECT(dimnames);
1820 ihaka 293 
	setAttrib(x, R_DimNamesSymbol, R_NilValue);
294 
	setAttrib(x, R_DimSymbol, newdims);
2506 maechler 295 
	if (dimnames != R_NilValue)
2715 pd 296 
	{
6994 pd 297 
	    if(!isNull(dnn))
298 
		setAttrib(newnames, R_NamesSymbol, newnamesnames);
2506 maechler 299 
	    setAttrib(x, R_DimNamesSymbol, newnames);
6994 pd 300 
	    UNPROTECT(2);
2715 pd 301 
	}
1858 ihaka 302 
	UNPROTECT(2);
1820 ihaka 303 
    }
304 
    UNPROTECT(1);
305 
    return x;
2 r 306 
}
307 
 
308 
SEXP do_drop(SEXP call, SEXP op, SEXP args, SEXP rho)
309 
{
1820 ihaka 310 
    SEXP x, xdims;
311 
    int i, n, shorten;
2 r 312 
 
1820 ihaka 313 
    checkArity(op, args);
314 
    x = CAR(args);
1839 ihaka 315 
    if ((xdims = getAttrib(x, R_DimSymbol)) != R_NilValue) {
1820 ihaka 316 
	n = LENGTH(xdims);
317 
	shorten = 0;
1839 ihaka 318 
	for (i = 0; i < n; i++)
319 
	    if (INTEGER(xdims)[i] == 1) shorten = 1;
320 
	if (shorten) {
321 
	    if (NAMED(x)) x = duplicate(x);
1820 ihaka 322 
	    x = DropDims(x);
2 r 323 
	}
1820 ihaka 324 
    }
325 
    return x;
2 r 326 
}
327 
 
1820 ihaka 328 
/* Length of Primitive Objects */
2 r 329 
 
330 
SEXP do_length(SEXP call, SEXP op, SEXP args, SEXP rho)
331 
{
1820 ihaka 332 
    SEXP ans;
28433 ripley 333 
    R_len_t len;
15783 rgentlem 334 
 
1820 ihaka 335 
    if (length(args) != 1)
5731 ripley 336 
	error("incorrect number of args to length");
15783 rgentlem 337 
 
16707 jmc 338 
    if( isObject(CAR(args)) && DispatchOrEval(call, op, "length", args,
17515 ripley 339 
					      rho, &ans, 0, 1))
15783 rgentlem 340 
      return(ans);
17515 ripley 341 
 
1820 ihaka 342 
    ans = allocVector(INTSXP, 1);
28433 ripley 343 
    len = length(CAR(args));
344 
    INTEGER(ans)[0] = (len < INT_MAX) ? len : NA_INTEGER;
1820 ihaka 345 
    return ans;
2 r 346 
}
347 
 
348 
 
349 
SEXP do_rowscols(SEXP call, SEXP op, SEXP args, SEXP rho)
350 
{
1820 ihaka 351 
    SEXP ans;
352 
    int i, j, nr, nc;
2 r 353 
 
1820 ihaka 354 
    if (length(args) != 1)
5731 ripley 355 
	error("incorrect number of args to row/col");
1820 ihaka 356 
    if (!isMatrix(CAR(args)))
5731 ripley 357 
	error("a matrix is required as arg to row/col");
2 r 358 
 
1820 ihaka 359 
    nr = nrows(CAR(args));
360 
    nc = ncols(CAR(args));
2 r 361 
 
1820 ihaka 362 
    ans = allocMatrix(INTSXP, nr, nc);
2 r 363 
 
1820 ihaka 364 
    switch (PRIMVAL(op)) {
365 
    case 1:
366 
	for (i = 0; i < nr; i++)
367 
	    for (j = 0; j < nc; j++)
368 
		INTEGER(ans)[i + j * nr] = i + 1;
369 
	break;
370 
    case 2:
371 
	for (i = 0; i < nr; i++)
372 
	    for (j = 0; j < nc; j++)
373 
		INTEGER(ans)[i + j * nr] = j + 1;
374 
	break;
375 
    }
376 
    return ans;
2 r 377 
}
378 
 
13942 bates 379 
static void matprod(double *x, int nrx, int ncx,
380 
		    double *y, int nry, int ncy, double *z)
381 
{
382 
    char *transa = "N", *transb = "N";
27297 ripley 383 
    int i,  j, k;
384 
    double one = 1.0, zero = 0.0, sum;
385 
    Rboolean have_na = FALSE;
386 
 
13996 duncan 387 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
27297 ripley 388 
	/* Don't trust the BLAS to handle NA/NaNs correctly: PR#4582
389 
	 * The test is only O(n) here
390 
	 */
391 
	for (i = 0; i < nrx*ncx; i++)
392 
	    if (ISNAN(x[i])) {have_na = TRUE; break;}
393 
	if (!have_na)
394 
	    for (i = 0; i < nry*ncy; i++)
395 
		if (ISNAN(y[i])) {have_na = TRUE; break;}
396 
	if (have_na) {
397 
	    for (i = 0; i < nrx; i++)
398 
		for (k = 0; k < ncy; k++) {
399 
		    sum = 0.0;
400 
		    for (j = 0; j < ncx; j++)
401 
			sum += x[i + j * nrx] * y[j + k * nry];
402 
		    z[i + k * nrx] = sum;
403 
		}
404 
	} else
405 
	    F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
406 
			    x, &nrx, y, &nry, &zero, z, &nrx);
407 
    } else /* zero-extent operations should return zeroes */
23175 ripley 408 
	for(i = 0; i < nrx*ncy; i++) z[i] = 0;
27297 ripley 409 
}
571 ihaka 410 
 
23236 ripley 411 
#ifdef HAVE_DOUBLE_COMPLEX
412 
/* ZGEMM - perform one of the matrix-matrix operations    */
23175 ripley 413 
/* C := alpha*op( A )*op( B ) + beta*C */
26467 maechler 414 
extern void
23175 ripley 415 
F77_NAME(zgemm)(const char *transa, const char *transb, const int *m,
416 
		const int *n, const int *k, const Rcomplex *alpha,
417 
		const Rcomplex *a, const int *lda,
418 
		const Rcomplex *b, const int *ldb,
419 
		const Rcomplex *beta, Rcomplex *c, const int *ldc);
23236 ripley 420 
#endif
23175 ripley 421 
 
6994 pd 422 
static void cmatprod(Rcomplex *x, int nrx, int ncx,
23175 ripley 423 
		     Rcomplex *y, int nry, int ncy, Rcomplex *z)
2 r 424 
{
29410 ripley 425 
#ifdef HAVE_DOUBLE_COMPLEX
23175 ripley 426 
    char *transa = "N", *transb = "N";
427 
    int i;
428 
    Rcomplex one, zero;
429 
 
430 
    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
431 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
432 
        F77_CALL(zgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
433 
			x, &nrx, y, &nry, &zero, z, &nrx);
434 
    } else { /* zero-extent operations should return zeroes */
435 
	for(i = 0; i < nrx*ncy; i++) z[i].r = z[i].i = 0;
436 
    }
437 
#else
1820 ihaka 438 
    int i, j, k;
439 
    double xij_r, xij_i, yjk_r, yjk_i, sum_i, sum_r;
2 r 440 
 
1839 ihaka 441 
    for (i = 0; i < nrx; i++)
442 
	for (k = 0; k < ncy; k++) {
443 
	    z[i + k * nrx].r = NA_REAL;
444 
	    z[i + k * nrx].i = NA_REAL;
1820 ihaka 445 
	    sum_r = 0.0;
446 
	    sum_i = 0.0;
1839 ihaka 447 
	    for (j = 0; j < ncx; j++) {
448 
		xij_r = x[i + j * nrx].r;
449 
		xij_i = x[i + j * nrx].i;
450 
		yjk_r = y[j + k * nry].r;
451 
		yjk_i = y[j + k * nry].i;
1820 ihaka 452 
		if (ISNAN(xij_r) || ISNAN(xij_i)
453 
		    || ISNAN(yjk_r) || ISNAN(yjk_i))
454 
		    goto next_ik;
455 
		sum_r += (xij_r * yjk_r - xij_i * yjk_i);
456 
		sum_i += (xij_r * yjk_i + xij_i * yjk_r);
457 
	    }
1839 ihaka 458 
	    z[i + k * nrx].r = sum_r;
459 
	    z[i + k * nrx].i = sum_i;
1820 ihaka 460 
	next_ik:
461 
	    ;
23175 ripley 462 
	}
1016 maechler 463 
#endif
2 r 464 
}
465 
 
18753 ripley 466 
static void symcrossprod(double *x, int nr, int nc, double *z)
467 
{
468 
    char *trans = "T", *uplo = "U";
469 
    double one = 1.0, zero = 0.0;
470 
    int i, j;
471 
    if (nr > 0 && nc > 0) {
472 
        F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc);
26467 maechler 473 
	for (i = 1; i < nc; i++)
18753 ripley 474 
	    for (j = 0; j < i; j++) z[i + nc *j] = z[j + nc * i];
26666 ripley 475 
    } else { /* zero-extent operations should return zeroes */
476 
	for(i = 0; i < nc*nc; i++) z[i] = 0;
18753 ripley 477 
    }
26666 ripley 478 
 
18753 ripley 479 
}
480 
 
1820 ihaka 481 
static void crossprod(double *x, int nrx, int ncx,
482 
		      double *y, int nry, int ncy, double *z)
2 r 483 
{
13942 bates 484 
    char *transa = "T", *transb = "N";
485 
    double one = 1.0, zero = 0.0;
13996 duncan 486 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
487 
        F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
23175 ripley 488 
			x, &nrx, y, &nry, &zero, z, &ncx);
13996 duncan 489 
    }
26666 ripley 490 
    else { /* zero-extent operations should return zeroes */
491 
	int i;
492 
	for(i = 0; i < ncx*ncy; i++) z[i] = 0;
493 
    }
2 r 494 
}
495 
 
6994 pd 496 
static void ccrossprod(Rcomplex *x, int nrx, int ncx,
497 
		       Rcomplex *y, int nry, int ncy, Rcomplex *z)
2 r 498 
{
23175 ripley 499 
    char *transa = "T", *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) {
504 
        F77_CALL(zgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
505 
			x, &nrx, y, &nry, &zero, z, &ncx);
506 
    }
26666 ripley 507 
    else { /* zero-extent operations should return zeroes */
508 
	int i;
509 
	for(i = 0; i < ncx*ncy; i++) z[i].r = z[i].i = 0;
510 
    }
2 r 511 
}
26467 maechler 512 
/* "%*%" (op = 0)  or  crossprod (op = 1) : */
2 r 513 
SEXP do_matprod(SEXP call, SEXP op, SEXP args, SEXP rho)
514 
{
1820 ihaka 515 
    int ldx, ldy, nrx, ncx, nry, ncy, mode;
18753 ripley 516 
    SEXP x = CAR(args), y = CADR(args), xdims, ydims, ans;
517 
    Rboolean sym;
2 r 518 
 
21875 jmc 519 
    if(R_has_methods(op)) {
520 
      SEXP value;
521 
      value = R_possible_dispatch(call, op, args, rho);
522 
      if(value) return value;
523 
    }
524 
 
18753 ripley 525 
    sym = isNull(y);
526 
    if (sym && (PRIMVAL(op) == 1)) y = x;
527 
    if ( !(isNumeric(x) || isComplex(x)) || !(isNumeric(y) || isComplex(y)) )
6206 rgentlem 528 
	errorcall(call, "requires numeric matrix/vector arguments");
2 r 529 
 
1820 ihaka 530 
    xdims = getAttrib(x, R_DimSymbol);
531 
    ydims = getAttrib(y, R_DimSymbol);
532 
    ldx = length(xdims);
533 
    ldy = length(ydims);
2 r 534 
 
26472 ihaka 535 
    if (ldx != 2 && ldy != 2) {		/* x and y non-matrices */
1839 ihaka 536 
	if (PRIMVAL(op) == 0) {
1820 ihaka 537 
	    nrx = 1;
538 
	    ncx = LENGTH(x);
2 r 539 
	}
1820 ihaka 540 
	else {
541 
	    nrx = LENGTH(x);
542 
	    ncx = 1;
2 r 543 
	}
1820 ihaka 544 
	nry = LENGTH(y);
545 
	ncy = 1;
546 
    }
26472 ihaka 547 
    else if (ldx != 2) {		/* x not a matrix */
1820 ihaka 548 
	nry = INTEGER(ydims)[0];
549 
	ncy = INTEGER(ydims)[1];
550 
	nrx = 0;
551 
	ncx = 0;
1839 ihaka 552 
	if (PRIMVAL(op) == 0) {
26472 ihaka 553 
	    if (LENGTH(x) == nry) {	/* x as row vector */
1820 ihaka 554 
		nrx = 1;
555 
		ncx = LENGTH(x);
2213 maechler 556 
	    }
26472 ihaka 557 
	    else if (nry == 1) {	/* x as col vector */
1820 ihaka 558 
		nrx = LENGTH(x);
559 
		ncx = 1;
560 
	    }
2 r 561 
	}
562 
	else {
26472 ihaka 563 
	    if (LENGTH(x) == nry) {	/* x is a row vector */
1820 ihaka 564 
		nrx = LENGTH(x);
565 
		ncx = 1;
566 
	    }
2 r 567 
	}
1820 ihaka 568 
    }
26472 ihaka 569 
    else if (ldy != 2) {		/* y not a matrix */
1820 ihaka 570 
	nrx = INTEGER(xdims)[0];
571 
	ncx = INTEGER(xdims)[1];
572 
	nry = 0;
573 
	ncy = 0;
1839 ihaka 574 
	if (PRIMVAL(op) == 0) {
26472 ihaka 575 
	    if (LENGTH(y) == ncx) {	/* y as col vector */
1820 ihaka 576 
		nry = LENGTH(y);
577 
		ncy = 1;
578 
	    }
26472 ihaka 579 
	    else if (ncx == 1) {	/* y as row vector */
1820 ihaka 580 
		nry = 1;
581 
		ncy = LENGTH(y);
582 
	    }
2 r 583 
	}
584 
	else {
26472 ihaka 585 
	    if (LENGTH(y) == nrx) {	/* y is a row vector */
1820 ihaka 586 
		nry = LENGTH(y);
587 
		ncy = 1;
588 
	    }
2 r 589 
	}
1820 ihaka 590 
    }
26472 ihaka 591 
    else {				/* x and y matrices */
1820 ihaka 592 
	nrx = INTEGER(xdims)[0];
593 
	ncx = INTEGER(xdims)[1];
594 
	nry = INTEGER(ydims)[0];
595 
	ncy = INTEGER(ydims)[1];
596 
    }
2 r 597 
 
1839 ihaka 598 
    if (PRIMVAL(op) == 0) {
599 
	if (ncx != nry)
5731 ripley 600 
	    errorcall(call, "non-conformable arguments");
1820 ihaka 601 
    }
602 
    else {
1839 ihaka 603 
	if (nrx != nry)
5731 ripley 604 
	    errorcall(call, "non-conformable arguments");
1820 ihaka 605 
    }
606 
 
1839 ihaka 607 
    if (isComplex(CAR(args)) || isComplex(CADR(args)))
1820 ihaka 608 
	mode = CPLXSXP;
609 
    else
610 
	mode = REALSXP;
10172 luke 611 
    SETCAR(args, coerceVector(CAR(args), mode));
612 
    SETCADR(args, coerceVector(CADR(args), mode));
1820 ihaka 613 
 
26472 ihaka 614 
    if (PRIMVAL(op) == 0) {		       	/* op == 0 : matprod() */
615 
 
1820 ihaka 616 
	PROTECT(ans = allocMatrix(mode, nrx, ncy));
1839 ihaka 617 
	if (mode == CPLXSXP)
1820 ihaka 618 
	    cmatprod(COMPLEX(CAR(args)), nrx, ncx,
619 
		     COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));
2 r 620 
	else
1820 ihaka 621 
	    matprod(REAL(CAR(args)), nrx, ncx,
622 
		    REAL(CADR(args)), nry, ncy, REAL(ans));
26472 ihaka 623 
 
1820 ihaka 624 
	PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));
625 
	PROTECT(ydims = getAttrib(CADR(args), R_DimNamesSymbol));
26472 ihaka 626 
 
1820 ihaka 627 
	if (xdims != R_NilValue || ydims != R_NilValue) {
6994 pd 628 
	    SEXP dimnames, dimnamesnames, dn;
629 
	    PROTECT(dimnames = allocVector(VECSXP, 2));
630 
	    PROTECT(dimnamesnames = allocVector(STRSXP, 2));
631 
	    if (xdims != R_NilValue) {
8628 pd 632 
		if (ldx == 2 || ncx ==1) {
633 
		    dn = getAttrib(xdims, R_NamesSymbol);
10172 luke 634 
		    SET_VECTOR_ELT(dimnames, 0, VECTOR_ELT(xdims, 0));
8628 pd 635 
		    if(!isNull(dn))
10172 luke 636 
			SET_STRING_ELT(dimnamesnames, 0, STRING_ELT(dn, 0));
8628 pd 637 
		}
6994 pd 638 
	    }
639 
	    if (ydims != R_NilValue) {
8628 pd 640 
		if (ldy == 2 ){
641 
		    dn = getAttrib(ydims, R_NamesSymbol);
10172 luke 642 
		    SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 1));
8628 pd 643 
		    if(!isNull(dn))
10172 luke 644 
			SET_STRING_ELT(dimnamesnames, 1, STRING_ELT(dn, 1));
8628 pd 645 
		} else if (nry == 1) {
646 
		    dn = getAttrib(ydims, R_NamesSymbol);
10172 luke 647 
		    SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 0));
8628 pd 648 
		    if(!isNull(dn))
10172 luke 649 
			SET_STRING_ELT(dimnamesnames, 1, STRING_ELT(dn, 0));
8628 pd 650 
		}
6994 pd 651 
	    }
652 
	    setAttrib(dimnames, R_NamesSymbol, dimnamesnames);
1858 ihaka 653 
	    setAttrib(ans, R_DimNamesSymbol, dimnames);
6994 pd 654 
	    UNPROTECT(2);
2 r 655 
	}
1820 ihaka 656 
    }
26472 ihaka 657 
 
658 
    else {					/* op == 1: crossprod() */
659 
 
1820 ihaka 660 
	PROTECT(ans = allocMatrix(mode, ncx, ncy));
1839 ihaka 661 
	if (mode == CPLXSXP)
23175 ripley 662 
	    if(sym)
663 
		ccrossprod(COMPLEX(CAR(args)), nrx, ncx,
664 
			   COMPLEX(CAR(args)), nry, ncy, COMPLEX(ans));
665 
	    else
666 
		ccrossprod(COMPLEX(CAR(args)), nrx, ncx,
667 
			   COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));
18753 ripley 668 
	else {
669 
	    if(sym)
670 
		symcrossprod(REAL(CAR(args)), nrx, ncx, REAL(ans));
671 
	    else
29387 ripley 672 
		crossprod(REAL(CAR(args)), nrx, ncx,
673 
			  REAL(CADR(args)), nry, ncy, REAL(ans));
18753 ripley 674 
	}
26472 ihaka 675 
 
1820 ihaka 676 
	PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));
26472 ihaka 677 
	if (sym)
18775 ripley 678 
	    PROTECT(ydims = xdims);
679 
	else
680 
	    PROTECT(ydims = getAttrib(CADR(args), R_DimNamesSymbol));
26472 ihaka 681 
 
1820 ihaka 682 
	if (xdims != R_NilValue || ydims != R_NilValue) {
7081 pd 683 
	    SEXP dimnames, dimnamesnames, dnx=R_NilValue, dny=R_NilValue;
26472 ihaka 684 
 
685 
	    /* allocate dimnames and dimnamesnames */
686 
 
6994 pd 687 
	    PROTECT(dimnames = allocVector(VECSXP, 2));
688 
	    PROTECT(dimnamesnames = allocVector(STRSXP, 2));
26472 ihaka 689 
 
690 
            /* There was a bug here.  The second element of a */
691 
            /* dimnames list was being accessed for a 1-d array. */
692 
            /* I have just excluded the use of dimnames in this */
693 
            /* case. - ihaka Sep 30, 2003. */
694 
 
6994 pd 695 
	    if (xdims != R_NilValue) {
26472 ihaka 696 
                if (ldx == 2) {
697 
		    dnx = getAttrib(xdims, R_NamesSymbol);
698 
		    SET_VECTOR_ELT(dimnames, 0, VECTOR_ELT(xdims, 1));
699 
		    if(!isNull(dnx))
700 
			SET_STRING_ELT(dimnamesnames, 0, STRING_ELT(dnx, 1));
701 
		}
6994 pd 702 
	    }
26472 ihaka 703 
 
6994 pd 704 
	    if (ydims != R_NilValue) {
26472 ihaka 705 
                if (ldy == 2) {
706 
		    dny = getAttrib(ydims, R_NamesSymbol);
707 
		    SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 1));
708 
		    if(!isNull(dny))
709 
			SET_STRING_ELT(dimnamesnames, 1, STRING_ELT(dny, 1));
710 
		}
6994 pd 711 
	    }
26472 ihaka 712 
 
713 
            /* We sometimes attach a dimnames attribute */
714 
            /* whose elements are all NULL ... */
715 
            /* Thus is ugly but causes no real damage. */
716 
 
7081 pd 717 
	    if (!isNull(dnx) || !isNull(dny))
6994 pd 718 
		setAttrib(dimnames, R_NamesSymbol, dimnamesnames);
1858 ihaka 719 
	    setAttrib(ans, R_DimNamesSymbol, dimnames);
6994 pd 720 
	    UNPROTECT(2);
2 r 721 
	}
1820 ihaka 722 
    }
723 
    UNPROTECT(3);
724 
    return ans;
2 r 725 
}
726 
 
727 
SEXP do_transpose(SEXP call, SEXP op, SEXP args, SEXP rho)
728 
{
9147 maechler 729 
    SEXP a, r, dims, dimnames, dimnamesnames=R_NilValue,
6994 pd 730 
	ndimnamesnames, rnames, cnames;
1839 ihaka 731 
    int i, len = 0, ncol=0, nrow=0;
2 r 732 
 
1820 ihaka 733 
    checkArity(op, args);
734 
    a = CAR(args);
2 r 735 
 
1839 ihaka 736 
    if (isVector(a)) {
1820 ihaka 737 
	dims = getAttrib(a, R_DimSymbol);
1858 ihaka 738 
	rnames = R_NilValue;
739 
	cnames = R_NilValue;
1820 ihaka 740 
	switch(length(dims)) {
741 
	case 0:
1858 ihaka 742 
	    nrow = len = length(a);
743 
	    ncol = 1;
744 
	    rnames = getAttrib(a, R_NamesSymbol);
745 
	    break;
1820 ihaka 746 
	case 1:
747 
	    nrow = len = length(a);
748 
	    ncol = 1;
1858 ihaka 749 
	    rnames = getAttrib(a, R_DimNamesSymbol);
750 
	    if (rnames != R_NilValue)
10172 luke 751 
		rnames = VECTOR_ELT(rnames, 0);
1820 ihaka 752 
	    break;
753 
	case 2:
754 
	    ncol = ncols(a);
755 
	    nrow = nrows(a);
756 
	    len = length(a);
1858 ihaka 757 
	    dimnames = getAttrib(a, R_DimNamesSymbol);
758 
	    if (dimnames != R_NilValue) {
10172 luke 759 
		rnames = VECTOR_ELT(dimnames, 0);
760 
		cnames = VECTOR_ELT(dimnames, 1);
6994 pd 761 
		dimnamesnames = getAttrib(dimnames, R_NamesSymbol);
1858 ihaka 762 
	    }
1820 ihaka 763 
	    break;
764 
	default:
765 
	    goto not_matrix;
2 r 766 
	}
1820 ihaka 767 
    }
1858 ihaka 768 
    else
769 
	goto not_matrix;
1820 ihaka 770 
    PROTECT(r = allocVector(TYPEOF(a), len));
771 
    switch (TYPEOF(a)) {
772 
    case LGLSXP:
773 
    case INTSXP:
774 
	for (i = 0; i < len; i++)
775 
	    INTEGER(r)[i] = INTEGER(a)[(i / ncol) + (i % ncol) * nrow];
776 
	break;
777 
    case REALSXP:
778 
	for (i = 0; i < len; i++)
779 
	    REAL(r)[i] = REAL(a)[(i / ncol) + (i % ncol) * nrow];
780 
	break;
781 
    case CPLXSXP:
782 
	for (i = 0; i < len; i++)
783 
	    COMPLEX(r)[i] = COMPLEX(a)[(i / ncol) + (i % ncol) * nrow];
784 
	break;
785 
    case STRSXP:
786 
	for (i = 0; i < len; i++)
10172 luke 787 
	    SET_STRING_ELT(r, i,
788 
			   STRING_ELT(a, (i / ncol) + (i % ncol) * nrow));
1820 ihaka 789 
	break;
1858 ihaka 790 
    case VECSXP:
791 
	for (i = 0; i < len; i++)
10172 luke 792 
	    SET_VECTOR_ELT(r, i,
793 
			   VECTOR_ELT(a, (i / ncol) + (i % ncol) * nrow));
1858 ihaka 794 
	break;
30691 murdoch 795 
    case RAWSXP:
796 
	for (i = 0; i < len; i++)
797 
	    RAW(r)[i] = RAW(a)[(i / ncol) + (i % ncol) * nrow];
798 
	break;
1858 ihaka 799 
    default:
800 
	goto not_matrix;
1820 ihaka 801 
    }
1858 ihaka 802 
    PROTECT(dims = allocVector(INTSXP, 2));
1820 ihaka 803 
    INTEGER(dims)[0] = ncol;
804 
    INTEGER(dims)[1] = nrow;
805 
    setAttrib(r, R_DimSymbol, dims);
1858 ihaka 806 
    UNPROTECT(1);
807 
    if(rnames != R_NilValue || cnames != R_NilValue) {
808 
	PROTECT(dimnames = allocVector(VECSXP, 2));
10172 luke 809 
	SET_VECTOR_ELT(dimnames, 0, cnames);
810 
	SET_VECTOR_ELT(dimnames, 1, rnames);
6994 pd 811 
	if(!isNull(dimnamesnames)) {
812 
	    PROTECT(ndimnamesnames = allocVector(VECSXP, 2));
10172 luke 813 
	    SET_STRING_ELT(ndimnamesnames, 1, STRING_ELT(dimnamesnames, 0));
814 
	    SET_STRING_ELT(ndimnamesnames, 0, STRING_ELT(dimnamesnames, 1));
6994 pd 815 
	    setAttrib(dimnames, R_NamesSymbol, ndimnamesnames);
816 
	    UNPROTECT(1);
817 
	}
1858 ihaka 818 
	setAttrib(r, R_DimNamesSymbol, dimnames);
2 r 819 
	UNPROTECT(1);
1820 ihaka 820 
    }
1858 ihaka 821 
    copyMostAttrib(a, r);
1820 ihaka 822 
    UNPROTECT(1);
823 
    return r;
824 
 not_matrix:
5731 ripley 825 
    errorcall(call, "argument is not a matrix");
1820 ihaka 826 
    return call;/* never used; just for -Wall */
2 r 827 
}
828 
 
13214 maechler 829 
/*
830 
 New version of aperm, using strides for speed.
831 
 Jonathan Rougier <J.C.Rougier@durham.ac.uk>
832 
 
833 
 v1.0 30.01.01
834 
 
835 
 M.Maechler : expanded	all ../include/Rdefines.h macros
836 
 */
837 
 
838 
/* this increments iip and sets j using strides */
839 
 
840 
#define CLICKJ						\
841 
    for (itmp=0; itmp<n; itmp++)			\
842 
	if (iip[itmp] == INTEGER(dimsr)[itmp]-1)	\
843 
	    iip[itmp] = 0;				\
844 
	else {						\
845 
	    iip[itmp]++;				\
846 
	    break;					\
847 
	}						\
848 
    for (j=0, itmp=0; itmp<n; itmp++)			\
849 
	j += iip[itmp] * stride[itmp];
17515 ripley 850 
 
13214 maechler 851 
/* aperm (a, perm, resize = TRUE) */
2 r 852 
SEXP do_aperm(SEXP call, SEXP op, SEXP args, SEXP rho)
853 
{
13214 maechler 854 
    SEXP a, perm, resize, r, dimsa, dimsr, dna;
855 
    int i, j, n, len, itmp;
856 
    int *pp, *iip, *stride;
857 
    char *vmax;
2 r 858 
 
1820 ihaka 859 
    checkArity(op, args);
2 r 860 
 
1820 ihaka 861 
    a = CAR(args);
13214 maechler 862 
    if (!isArray(a))
863 
	errorcall(call,"invalid first argument, must be an array");
864 
 
1820 ihaka 865 
    PROTECT(dimsa = getAttrib(a, R_DimSymbol));
13214 maechler 866 
    n = LENGTH(dimsa);
2 r 867 
 
13214 maechler 868 
    /* check the permutation */
869 
 
1820 ihaka 870 
    PROTECT(perm = coerceVector(CADR(args), INTSXP));
13214 maechler 871 
    vmax = vmaxget();
872 
    pp = (int *) R_alloc(n, sizeof(int));
873 
    if (length(perm) == 0) {
874 
	for (i=0; i<n; i++)
875 
	    pp[i] = n-1-i;
876 
    } else if (length(perm) == n) {
877 
	for (i=0; i<n; i++)
878 
	    pp[i] = INTEGER(perm)[i] - 1; /* no offset! */
879 
    } else
880 
	errorcall(call, "`perm' is of wrong length");
2 r 881 
 
13214 maechler 882 
    iip = (int *) R_alloc(n, sizeof(int));
883 
    for (i=0; i<n; iip[i++] = 0);
884 
    for (i=0; i<n; i++)
885 
	if (pp[i] >= 0 && pp[i] < n)
886 
	    iip[pp[i]]++;
887 
	else
888 
	    errorcall(call, "value of out range in `perm'");
889 
    for (i=0; i<n; i++)
890 
	if (iip[i]==0)
891 
	    errorcall(call, "invalid permutation (`perm')");
2 r 892 
 
13214 maechler 893 
    /* create the stride object and permute */
2 r 894 
 
13214 maechler 895 
    stride = (int *) R_alloc(n, sizeof(int));
896 
 
897 
    for (iip[0] = 1, i = 1; i<n; i++)
898 
	iip[i] = iip[i-1] * INTEGER(dimsa)[i-1];
899 
 
900 
    for (i=0; i<n; i++)
901 
	stride[i] = iip[pp[i]];
902 
 
903 
    /* also need to have the dimensions of r */
904 
 
905 
    PROTECT(dimsr = allocVector(INTSXP,n));
906 
    for (i=0; i<n; i++)
907 
	INTEGER(dimsr)[i] = INTEGER(dimsa)[pp[i]];
908 
 
909 
    /* and away we go! iip will hold the incrementer */
910 
 
911 
    len = LENGTH(a);
912 
    len = length(a);
1820 ihaka 913 
    PROTECT(r = allocVector(TYPEOF(a), len));
2 r 914 
 
13214 maechler 915 
    for (i=0; i<n; iip[i++] = 0);
916 
 
1820 ihaka 917 
    switch (TYPEOF(a)) {
13214 maechler 918 
 
1820 ihaka 919 
    case INTSXP:
920 
    case LGLSXP:
13214 maechler 921 
	for (j=0, i=0; i<len; i++) {
922 
	    INTEGER(r)[i] = INTEGER(a)[j];
923 
	    CLICKJ;
2 r 924 
	}
1820 ihaka 925 
	break;
13214 maechler 926 
 
1820 ihaka 927 
    case REALSXP:
13214 maechler 928 
	for (j=0, i=0; i<len; i++) {
929 
	    REAL(r)[i] = REAL(a)[j];
930 
	    CLICKJ;
1820 ihaka 931 
	}
932 
	break;
13214 maechler 933 
 
1820 ihaka 934 
    case CPLXSXP:
13214 maechler 935 
	for (j=0, i=0; i<len; i++) {
936 
	    COMPLEX(r)[i].r = COMPLEX(a)[j].r;
937 
	    COMPLEX(r)[i].i = COMPLEX(a)[j].i;
938 
	    CLICKJ;
1820 ihaka 939 
	}
940 
	break;
13214 maechler 941 
 
1820 ihaka 942 
    case STRSXP:
13214 maechler 943 
	for (j=0, i=0; i<len; i++) {
944 
	    SET_STRING_ELT(r, i, STRING_ELT(a, j));
945 
	    CLICKJ;
1820 ihaka 946 
	}
947 
	break;
13214 maechler 948 
 
1839 ihaka 949 
    case VECSXP:
13214 maechler 950 
	for (j=0, i=0; i<len; i++) {
951 
	    SET_VECTOR_ELT(r, i, VECTOR_ELT(a, j));
952 
	    CLICKJ;
1839 ihaka 953 
	}
12976 pd 954 
	break;
13214 maechler 955 
 
30691 murdoch 956 
    case RAWSXP:
957 
	for (j=0, i=0; i<len; i++) {
958 
	    RAW(r)[i] = RAW(a)[j];
959 
	    CLICKJ;
960 
	}
961 
	break;
962 
 
1820 ihaka 963 
    default:
13214 maechler 964 
	errorcall(call, "unsupported type of array");
1820 ihaka 965 
    }
17515 ripley 966 
 
13214 maechler 967 
    /* handle the resize */
968 
    PROTECT(resize = coerceVector(CADDR(args), INTSXP));
969 
    if (LOGICAL(resize)[0])
1820 ihaka 970 
	setAttrib(r, R_DimSymbol, dimsr);
971 
    else
972 
	setAttrib(r, R_DimSymbol, dimsa);
13214 maechler 973 
 
974 
    /* and handle the dimnames */
975 
 
976 
    PROTECT(dna = getAttrib(a, R_DimNamesSymbol));
977 
 
978 
    if (LOGICAL(resize)[0] && dna != R_NilValue) {
979 
 
980 
	SEXP dnna, dnr, dnnr;
981 
 
982 
	PROTECT(dnna = getAttrib(dna, R_NamesSymbol));
983 
	PROTECT(dnnr = allocVector(STRSXP,n));
984 
	PROTECT(dnr  = allocVector(VECSXP,n));
985 
 
986 
	for (i=0; i<n; i++) {
987 
	    SET_VECTOR_ELT(dnr, i, VECTOR_ELT(dna, pp[i]));
988 
	    if (dnna != R_NilValue)
989 
		SET_STRING_ELT(dnnr, i, STRING_ELT(dnna, pp[i]));
990 
	}
991 
 
992 
	if (dnna != R_NilValue)
993 
	    setAttrib(dnr, R_NamesSymbol, dnnr);
994 
	setAttrib(r, R_DimNamesSymbol, dnr);
995 
	UNPROTECT(3); /* dnna, dnr, dnnr */
996 
    }
997 
    /* free temporary memory */
998 
    vmaxset(vmax);
17515 ripley 999 
 
13214 maechler 1000 
    UNPROTECT(6); /* dimsa, perm, r, dimsr, resize, dna */
1820 ihaka 1001 
    return r;
2 r 1002 
}
17515 ripley 1003 
 
1004 
/* colSums(x, n, p, na.rm) and friends */
1005 
SEXP do_colsum(SEXP call, SEXP op, SEXP args, SEXP rho)
1006 
{
1007 
    SEXP x, ans = R_NilValue;
1008 
    int OP, n, p, cnt = 0, i, j, type;
17521 ripley 1009 
    Rboolean NaRm, keepNA;
17515 ripley 1010 
    int *ix;
1011 
    double *rx, sum = 0.0;
1012 
 
1013 
    checkArity(op, args);
1014 
    x = CAR(args); args = CDR(args);
1015 
    n = asInteger(CAR(args)); args = CDR(args);
1016 
    p = asInteger(CAR(args)); args = CDR(args);
1017 
    NaRm = asLogical(CAR(args));
17521 ripley 1018 
    if (n == NA_INTEGER || n <= 0)
17515 ripley 1019 
	errorcall(call, "invalid value of n");
17521 ripley 1020 
    if (p == NA_INTEGER || p <= 0)
17515 ripley 1021 
	errorcall(call, "invalid value of p");
17521 ripley 1022 
    if (NaRm == NA_LOGICAL) errorcall(call, "invalid value of na.rm");
1023 
    keepNA = !NaRm;
17515 ripley 1024 
 
1025 
    OP = PRIMVAL(op);
1026 
    switch (type = TYPEOF(x)) {
1027 
    case LGLSXP: break;
1028 
    case INTSXP: break;
1029 
    case REALSXP: break;
1030 
    default:
1031 
	errorcall(call, "`x' must be numeric");
1032 
    }
1033 
 
17521 ripley 1034 
    if (OP == 0 || OP == 1) { /* columns */
1035 
	cnt = n;
17515 ripley 1036 
	PROTECT(ans = allocVector(REALSXP, p));
1037 
	for (j = 0; j < p; j++) {
1038 
	    switch (type) {
1039 
	    case REALSXP:
17521 ripley 1040 
		rx = REAL(x) + n*j;
1041 
		if (keepNA)
1042 
		    for (sum = 0., i = 0; i < n; i++) sum += *rx++;
1043 
		else {
1044 
		    for (cnt = 0, sum = 0., i = 0; i < n; i++, rx++)
1045 
			if (!ISNAN(*rx)) {cnt++; sum += *rx;}
1046 
			else if (keepNA) {sum = NA_REAL; break;}
1047 
		}
17515 ripley 1048 
		break;
1049 
	    case INTSXP:
17521 ripley 1050 
		ix = INTEGER(x) + n*j;
1051 
		for (cnt = 0, sum = 0., i = 0; i < n; i++, ix++)
1052 
		    if (*ix != NA_INTEGER) {cnt++; sum += *ix;}
1053 
		    else if (keepNA) {sum = NA_REAL; break;}
17515 ripley 1054 
		break;
1055 
	    case LGLSXP:
17521 ripley 1056 
		ix = LOGICAL(x) + n*j;
1057 
		for (cnt = 0, sum = 0., i = 0; i < n; i++, ix++)
1058 
		    if (*ix != NA_LOGICAL) {cnt++; sum += *ix;}
1059 
		    else if (keepNA) {sum = NA_REAL; break;}
17515 ripley 1060 
		break;
1061 
	    }
17521 ripley 1062 
	    if (OP == 1) {
1063 
		if (cnt > 0) sum /= cnt; else sum = NA_REAL;
17515 ripley 1064 
	    }
1065 
	    REAL(ans)[j] = sum;
1066 
	}
1067 
    }
1068 
 
17521 ripley 1069 
    if (OP == 2 || OP == 3) { /* rows */
1070 
	cnt = p;
17515 ripley 1071 
	PROTECT(ans = allocVector(REALSXP, n));
17521 ripley 1072 
 
1073 
	/* reverse summation order to improve cache hits */
1074 
	if (type == REALSXP) {
1075 
	    double *rans = REAL(ans), *ra = rans, *cnt = NULL, *c;
1076 
	    rx = REAL(x);
1077 
	    if (!keepNA && OP == 3) cnt = Calloc(n, double);
1078 
	    for (ra = rans, i = 0; i < n; i++) *ra++ = 0.0;
1079 
	    for (j = 0; j < p; j++) {
1080 
		ra = rans;
1081 
		if (keepNA)
1082 
		    for (i = 0; i < n; i++) *ra++ += *rx++;
1083 
		else
26467 maechler 1084 
		    for (c = cnt, i = 0; i < n; i++, ra++, rx++, c++)
17521 ripley 1085 
			if (!ISNAN(*rx)) {
1086 
			    *ra += *rx;
1087 
			    if (OP == 3) (*c)++;
1088 
			}
1089 
	    }
1090 
	    if (OP == 3) {
1091 
		if (keepNA)
26467 maechler 1092 
		    for (ra = rans, i = 0; i < n; i++)
17521 ripley 1093 
			*ra++ /= p;
1094 
		else {
26467 maechler 1095 
		    for (ra = rans, c = cnt, i = 0; i < n; i++, c++)
17521 ripley 1096 
			if (*c > 0) *ra++ /= *c; else *ra++ = NA_REAL;
1097 
		    Free(cnt);
1098 
		}
1099 
	    }
1100 
	    UNPROTECT(1);
1101 
	    return ans;
1102 
	}
1103 
 
17515 ripley 1104 
	for (i = 0; i < n; i++) {
1105 
	    switch (type) {
29410 ripley 1106 
	    case REALSXP: /* this cannot be reached */
17521 ripley 1107 
		rx = REAL(x) + i;
1108 
		if (keepNA)
1109 
		    for (sum = 0., j = 0; j < p; j++, rx += n) sum += *rx;
1110 
		else {
1111 
		    for (cnt = 0, sum = 0., j = 0; j < p; j++, rx += n)
1112 
			if (!ISNAN(*rx)) {cnt++; sum += *rx;}
1113 
			else if (keepNA) {sum = NA_REAL; break;}
1114 
		}
17515 ripley 1115 
		break;
1116 
	    case INTSXP:
17521 ripley 1117 
		ix = INTEGER(x) + i;
1118 
		for (cnt = 0, sum = 0., j = 0; j < p; j++, ix += n)
1119 
		    if (*ix != NA_INTEGER) {cnt++; sum += *ix;}
1120 
		    else if (keepNA) {sum = NA_REAL; break;}
17515 ripley 1121 
		break;
1122 
	    case LGLSXP:
17521 ripley 1123 
		ix = LOGICAL(x) + i;
1124 
		for (cnt = 0, sum = 0., j = 0; j < p; j++, ix += n)
1125 
		    if (*ix != NA_LOGICAL) {cnt++; sum += *ix;}
1126 
		    else if (keepNA) {sum = NA_REAL; break;}
17515 ripley 1127 
		break;
1128 
	    }
17521 ripley 1129 
	    if (OP == 3) {
1130 
		if (cnt > 0) sum /= cnt; else sum = NA_REAL;
17515 ripley 1131 
	    }
1132 
	    REAL(ans)[i] = sum;
1133 
	}
1134 
    }
1135 
 
1136 
    UNPROTECT(1);
1137 
    return ans;
1138 
}