The R Project SVN R

Rev

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