The R Project SVN R

Rev

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