The R Project SVN R

Rev

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