The R Project SVN R

Rev

Rev 5187 | Rev 5731 | 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
4 
 *
5 
 *  This program is free software; you can redistribute it and/or modify
6 
 *  it under the terms of the GNU General Public License as published by
7 
 *  the Free Software Foundation; either version 2 of the License, or
8 
 *  (at your option) any later version.
9 
 *
10 
 *  This program is distributed in the hope that it will be useful,
11 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
12 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 
 *  GNU General Public License for more details.
14 
 *
15 
 *  You should have received a copy of the GNU General Public License
16 
 *  along with this program; if not, write to the Free Software
5458 ripley 17 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
2 r 18 
 */
19 
 
5187 hornik 20 
#ifdef HAVE_CONFIG_H
21 
#include <Rconfig.h>
22 
#endif
23 
 
2 r 24 
#include "Defn.h"
25 
#include "Mathlib.h"
26 
 
1839 ihaka 27 
/* "GetRowNames" and "GetColNames" are utility routines which */
28 
/* locate and return the row names and column names from the */
29 
/* dimnames attribute of a matrix.  They are useful because */
30 
/* old versions of R used pair-based lists for dimnames */
31 
/* whereas recent versions use vector bassed lists */
32 
 
1881 ihaka 33 
/* FIXME : This is nonsense.  When the "dimnames" attribute is */
34 
/* grabbed off an array it is always adjusted to be a vector. */
35 
 
1839 ihaka 36 
SEXP GetRowNames(SEXP dimnames)
37 
{
38 
    if (TYPEOF(dimnames) == VECSXP)
39 
	return VECTOR(dimnames)[0];
40 
    else if (TYPEOF(dimnames) == LISTSXP)
41 
	return CAR(dimnames);
42 
    else
1858 ihaka 43 
	return R_NilValue;
2213 maechler 44 
}
1839 ihaka 45 
 
46 
SEXP GetColNames(SEXP dimnames)
47 
{
48 
    if (TYPEOF(dimnames) == VECSXP)
49 
	return VECTOR(dimnames)[1];
50 
    else if (TYPEOF(dimnames) == LISTSXP)
51 
	return CADR(dimnames);
52 
    else
2213 maechler 53 
	return R_NilValue;
1839 ihaka 54 
}
55 
 
2 r 56 
SEXP do_matrix(SEXP call, SEXP op, SEXP args, SEXP rho)
57 
{
1820 ihaka 58 
    SEXP vals, snr, snc;
59 
    int nr, nc, byrow, lendat;
2 r 60 
 
1820 ihaka 61 
    checkArity(op, args);
62 
    vals = CAR(args);
63 
    snr = CADR(args);
64 
    snc = CADDR(args);
65 
    byrow = asInteger(CADR(CDDR(args)));
2 r 66 
 
1820 ihaka 67 
    if (isVector(vals) || isList(vals)) {
1839 ihaka 68 
	if (length(vals) < 0)
1820 ihaka 69 
	    errorcall(call, "argument has length zero\n");
1881 ihaka 70 
    }
71 
    else errorcall(call, "invalid matrix element type\n");
2 r 72 
 
1820 ihaka 73 
    if (!isNumeric(snr) || !isNumeric(snc))
74 
	error("non-numeric matrix extent\n");
2 r 75 
 
1820 ihaka 76 
    lendat = length(vals);
77 
    nr = asInteger(snr);
78 
    nc = asInteger(snc);
2 r 79 
 
1839 ihaka 80 
    if (lendat > 1 && (nr * nc) % lendat != 0) {
81 
	if (((lendat > nr) && (lendat / nr) * nr != lendat) ||
82 
	    ((lendat < nr) && (nr / lendat) * lendat != nr))
1820 ihaka 83 
	    warning("Replacement length not a multiple of the elements to replace in matrix(...) \n");
1839 ihaka 84 
	else if (((lendat > nc) && (lendat / nc) * nc != lendat) ||
85 
		 ((lendat < nc) && (nc / lendat) * lendat != nc))
1820 ihaka 86 
	    warning("Replacement length not a multiple of the elements to replace in matrix(...) \n");
2213 maechler 87 
    }
1839 ihaka 88 
	else if ((lendat > 1) && (nr * nc == 0)){
1819 thomas 89 
	  warning("Replacement length not a multiple of the elements to replace in matrix(...) \n");
2 r 90 
	}
1839 ihaka 91 
	else if (lendat == 0 && nr * nc > 0){
1819 thomas 92 
	  error("No data to replace in matrix(...)\n");
93 
	}
2 r 94 
 
1820 ihaka 95 
    PROTECT(snr = allocMatrix(TYPEOF(vals), nr, nc));
1839 ihaka 96 
    if (isVector(vals))
1820 ihaka 97 
	copyMatrix(snr, vals, byrow);
98 
    else
99 
	copyListMatrix(snr, vals, byrow);
100 
    UNPROTECT(1);
101 
    return snr;
2 r 102 
}
103 
 
104 
 
105 
SEXP allocMatrix(SEXPTYPE mode, int nrow, int ncol)
106 
{
1820 ihaka 107 
    SEXP s, t;
108 
    int n;
2 r 109 
 
1820 ihaka 110 
    if (nrow < 0 || ncol < 0)
111 
	error("negative extents to matrix\n");
112 
    n = nrow * ncol;
113 
    PROTECT(s = allocVector(mode, n));
114 
    PROTECT(t = allocVector(INTSXP, 2));
115 
    INTEGER(t)[0] = nrow;
116 
    INTEGER(t)[1] = ncol;
117 
    setAttrib(s, R_DimSymbol, t);
118 
    UNPROTECT(2);
119 
    return s;
2 r 120 
}
121 
 
1130 maechler 122 
 
2 r 123 
SEXP allocArray(SEXPTYPE mode, SEXP dims)
124 
{
1820 ihaka 125 
    SEXP array;
126 
    int i, n;
2 r 127 
 
1820 ihaka 128 
    n = 1;
129 
    for (i = 0; i < LENGTH(dims); i++)
130 
	n = n * INTEGER(dims)[i];
2 r 131 
 
1820 ihaka 132 
    PROTECT(dims = duplicate(dims));
133 
    PROTECT(array = allocVector(mode, n));
134 
    setAttrib(array, R_DimSymbol, dims);
135 
    UNPROTECT(2);
136 
    return array;
2 r 137 
}
138 
 
1820 ihaka 139 
/* DropDims strips away redundant dimensioning information. */
140 
/* If there is an appropriate dimnames attribute the correct */
141 
/* element is extracted and attached to the vector as a names */
142 
/* attribute.  Note that this function mutates x. */
143 
/* Duplication should occur before this is called. */
2 r 144 
 
145 
SEXP DropDims(SEXP x)
146 
{
2506 maechler 147 
    SEXP q, dims, dimnames, newnames = R_NilValue;
1839 ihaka 148 
    int i, n, ndims;
2 r 149 
 
1820 ihaka 150 
    PROTECT(x);
151 
    dims = getAttrib(x, R_DimSymbol);
152 
    dimnames = getAttrib(x, R_DimNamesSymbol);
2 r 153 
 
1820 ihaka 154 
    /* Check that dropping will actually do something. */
155 
    /* (1) Check that there is a "dim" attribute. */
2 r 156 
 
1839 ihaka 157 
    if (dims == R_NilValue) {
1820 ihaka 158 
	UNPROTECT(1);
159 
	return x;
160 
    }
1839 ihaka 161 
    ndims = LENGTH(dims);
2 r 162 
 
1839 ihaka 163 
    /* (2) Check whether there are redundant extents */
1820 ihaka 164 
    n = 0;
1839 ihaka 165 
    for (i = 0; i < ndims; i++)
166 
	if (INTEGER(dims)[i] != 1) n++;
167 
    if (n == ndims) {
1820 ihaka 168 
	UNPROTECT(1);
169 
	return x;
170 
    }
2 r 171 
 
1839 ihaka 172 
    if (n <= 1) {
173 
	/* We have reduced to a vector result. */
1820 ihaka 174 
	if (dimnames != R_NilValue) {
175 
	    n = length(dims);
1839 ihaka 176 
	    if (TYPEOF(dimnames) == VECSXP) {
177 
		for (i = 0; i < n; i++) {
178 
		    if (INTEGER(dims)[i] != 1) {
179 
			newnames = VECTOR(dimnames)[i];
180 
			break;
181 
		    }
2 r 182 
		}
1820 ihaka 183 
	    }
1839 ihaka 184 
	    else {
185 
		q = dimnames;
186 
		for (i = 0; i < n; i++) {
187 
		    if (INTEGER(dims)[i] != 1) {
188 
			newnames = CAR(q);
189 
			break;
190 
		    }
191 
		    q = CDR(q);
1820 ihaka 192 
		}
193 
	    }
2 r 194 
	}
1820 ihaka 195 
	PROTECT(newnames);
196 
	setAttrib(x, R_DimNamesSymbol, R_NilValue);
197 
	setAttrib(x, R_DimSymbol, R_NilValue);
198 
	setAttrib(x, R_NamesSymbol, newnames);
199 
	UNPROTECT(1);
200 
    }
1839 ihaka 201 
    else {
202 
	/* We have a lower dimensional array. */
2506 maechler 203 
	SEXP newdims;
1820 ihaka 204 
	PROTECT(newdims = allocVector(INTSXP, n));
1858 ihaka 205 
	for (i = 0, n = 0; i < ndims; i++)
1839 ihaka 206 
	    if (INTEGER(dims)[i] != 1)
1820 ihaka 207 
		INTEGER(newdims)[n++] = INTEGER(dims)[i];
1858 ihaka 208 
	if (dimnames != R_NilValue) {
209 
	    int havenames = 0;
210 
	    for (i = 0; i < ndims; i++)
211 
		if (INTEGER(dims)[i] != 1 && VECTOR(dimnames)[i] != R_NilValue)
212 
		    havenames = 1;
213 
	    if (havenames) {
2715 pd 214 
		PROTECT(newnames = allocVector(VECSXP, n));
1858 ihaka 215 
		for (i = 0, n= 0; i < ndims; i++) {
216 
		    if (INTEGER(dims)[i] != 1)
2506 maechler 217 
			VECTOR(newnames)[n++] = VECTOR(dimnames)[i];
1858 ihaka 218 
		}
1820 ihaka 219 
	    }
1858 ihaka 220 
	    else dimnames = R_NilValue;
1839 ihaka 221 
	}
1858 ihaka 222 
	PROTECT(dimnames);
1820 ihaka 223 
	setAttrib(x, R_DimNamesSymbol, R_NilValue);
224 
	setAttrib(x, R_DimSymbol, newdims);
2506 maechler 225 
	if (dimnames != R_NilValue)
2715 pd 226 
	{
2506 maechler 227 
	    setAttrib(x, R_DimNamesSymbol, newnames);
2715 pd 228 
	    UNPROTECT(1);
229 
	}
1858 ihaka 230 
	UNPROTECT(2);
1820 ihaka 231 
    }
232 
    UNPROTECT(1);
233 
    return x;
2 r 234 
}
235 
 
236 
SEXP do_drop(SEXP call, SEXP op, SEXP args, SEXP rho)
237 
{
1820 ihaka 238 
    SEXP x, xdims;
239 
    int i, n, shorten;
2 r 240 
 
1820 ihaka 241 
    checkArity(op, args);
242 
    x = CAR(args);
1839 ihaka 243 
    if ((xdims = getAttrib(x, R_DimSymbol)) != R_NilValue) {
1820 ihaka 244 
	n = LENGTH(xdims);
245 
	shorten = 0;
1839 ihaka 246 
	for (i = 0; i < n; i++)
247 
	    if (INTEGER(xdims)[i] == 1) shorten = 1;
248 
	if (shorten) {
249 
	    if (NAMED(x)) x = duplicate(x);
1820 ihaka 250 
	    x = DropDims(x);
2 r 251 
	}
1820 ihaka 252 
    }
253 
    return x;
2 r 254 
}
255 
 
1820 ihaka 256 
/* Length of Primitive Objects */
2 r 257 
 
258 
SEXP do_length(SEXP call, SEXP op, SEXP args, SEXP rho)
259 
{
1820 ihaka 260 
    SEXP ans;
261 
    if (length(args) != 1)
262 
	error("incorrect number of args to length\n");
263 
    ans = allocVector(INTSXP, 1);
264 
    INTEGER(ans)[0] = length(CAR(args));
265 
    return ans;
2 r 266 
}
267 
 
268 
 
269 
SEXP do_rowscols(SEXP call, SEXP op, SEXP args, SEXP rho)
270 
{
1820 ihaka 271 
    SEXP ans;
272 
    int i, j, nr, nc;
2 r 273 
 
1820 ihaka 274 
    if (length(args) != 1)
275 
	error("incorrect number of args to row/col\n");
276 
    if (!isMatrix(CAR(args)))
277 
	error("a matrix is required as arg to row/col\n");
2 r 278 
 
1820 ihaka 279 
    nr = nrows(CAR(args));
280 
    nc = ncols(CAR(args));
2 r 281 
 
1820 ihaka 282 
    ans = allocMatrix(INTSXP, nr, nc);
2 r 283 
 
1820 ihaka 284 
    switch (PRIMVAL(op)) {
285 
    case 1:
286 
	for (i = 0; i < nr; i++)
287 
	    for (j = 0; j < nc; j++)
288 
		INTEGER(ans)[i + j * nr] = i + 1;
289 
	break;
290 
    case 2:
291 
	for (i = 0; i < nr; i++)
292 
	    for (j = 0; j < nc; j++)
293 
		INTEGER(ans)[i + j * nr] = j + 1;
294 
	break;
295 
    }
296 
    return ans;
2 r 297 
}
298 
 
1820 ihaka 299 
/* FIXME - What about non-IEEE overflow ??? */
300 
/* Does it really matter? */
571 ihaka 301 
 
1820 ihaka 302 
static void matprod(double *x, int nrx, int ncx,
303 
		    double *y, int nry, int ncy, double *z)
2 r 304 
{
1820 ihaka 305 
    int i, j, k;
306 
    double xij, yjk, sum;
2 r 307 
 
1820 ihaka 308 
    for (i = 0; i < nrx; i++)
309 
	for (k = 0; k < ncy; k++) {
310 
	    z[i + k * nrx] = NA_REAL;
311 
	    sum = 0.0;
312 
	    for (j = 0; j < ncx; j++) {
313 
		xij = x[i + j * nrx];
314 
		yjk = y[j + k * nry];
571 ihaka 315 
#ifndef IEEE_754
1820 ihaka 316 
		if (ISNAN(xij) || ISNAN(yjk))
317 
		    goto next_ik;
571 ihaka 318 
#endif
1820 ihaka 319 
		sum += xij * yjk;
320 
	    }
321 
	    z[i + k * nrx] = sum;
1016 maechler 322 
#ifndef IEEE_754
1820 ihaka 323 
	next_ik:
324 
	    ;
1016 maechler 325 
#endif
1820 ihaka 326 
	}
2 r 327 
}
328 
 
329 
static void cmatprod(complex *x, int nrx, int ncx,
330 
		complex *y, int nry, int ncy, complex *z)
331 
{
1820 ihaka 332 
    int i, j, k;
333 
    double xij_r, xij_i, yjk_r, yjk_i, sum_i, sum_r;
2 r 334 
 
1839 ihaka 335 
    for (i = 0; i < nrx; i++)
336 
	for (k = 0; k < ncy; k++) {
337 
	    z[i + k * nrx].r = NA_REAL;
338 
	    z[i + k * nrx].i = NA_REAL;
1820 ihaka 339 
	    sum_r = 0.0;
340 
	    sum_i = 0.0;
1839 ihaka 341 
	    for (j = 0; j < ncx; j++) {
342 
		xij_r = x[i + j * nrx].r;
343 
		xij_i = x[i + j * nrx].i;
344 
		yjk_r = y[j + k * nry].r;
345 
		yjk_i = y[j + k * nry].i;
571 ihaka 346 
#ifndef IEEE_754
1820 ihaka 347 
		if (ISNAN(xij_r) || ISNAN(xij_i)
348 
		    || ISNAN(yjk_r) || ISNAN(yjk_i))
349 
		    goto next_ik;
571 ihaka 350 
#endif
1820 ihaka 351 
		sum_r += (xij_r * yjk_r - xij_i * yjk_i);
352 
		sum_i += (xij_r * yjk_i + xij_i * yjk_r);
353 
	    }
1839 ihaka 354 
	    z[i + k * nrx].r = sum_r;
355 
	    z[i + k * nrx].i = sum_i;
1016 maechler 356 
#ifndef IEEE_754
1820 ihaka 357 
	next_ik:
358 
	    ;
1016 maechler 359 
#endif
1820 ihaka 360 
	}
2 r 361 
}
362 
 
1820 ihaka 363 
static void crossprod(double *x, int nrx, int ncx,
364 
		      double *y, int nry, int ncy, double *z)
2 r 365 
{
1820 ihaka 366 
    int i, j, k;
367 
    double xji, yjk, sum;
2 r 368 
 
1820 ihaka 369 
    for (i = 0; i < ncx; i++)
370 
	for (k = 0; k < ncy; k++) {
371 
	    z[i + k * ncx] = NA_REAL;
372 
	    sum = 0.0;
373 
	    for (j = 0; j < nrx; j++) {
374 
		xji = x[j + i * nrx];
375 
		yjk = y[j + k * nry];
571 ihaka 376 
#ifndef IEEE_754
1820 ihaka 377 
		if (ISNAN(xji) || ISNAN(yjk))
378 
		    goto next_ik;
571 ihaka 379 
#endif
1820 ihaka 380 
		sum += xji * yjk;
381 
	    }
382 
	    z[i + k * ncx] = sum;
1016 maechler 383 
#ifndef IEEE_754
1820 ihaka 384 
	next_ik:
385 
	    ;
1016 maechler 386 
#endif
1820 ihaka 387 
	}
2 r 388 
}
389 
 
1820 ihaka 390 
static void ccrossprod(complex *x, int nrx, int ncx,
391 
		       complex *y, int nry, int ncy, complex *z)
2 r 392 
{
1820 ihaka 393 
    int i, j, k;
394 
    double xji_r, xji_i, yjk_r, yjk_i, sum_r, sum_i;
2 r 395 
 
1820 ihaka 396 
    for (i = 0; i < ncx; i++)
397 
	for (k = 0; k < ncy; k++) {
398 
	    z[i + k * ncx].r = NA_REAL;
399 
	    z[i + k * ncx].i = NA_REAL;
400 
	    sum_r = 0.0;
401 
	    sum_i = 0.0;
402 
	    for (j = 0; j < nrx; j++) {
403 
		xji_r = x[j + i * nrx].r;
404 
		xji_i = x[j + i * nrx].i;
405 
		yjk_r = y[j + k * nry].r;
406 
		yjk_i = y[j + k * nry].i;
571 ihaka 407 
#ifndef IEEE_754
1820 ihaka 408 
		if (ISNAN(xji_r) || ISNAN(xji_i)
409 
		    || ISNAN(yjk_r) || ISNAN(yjk_i))
410 
		    goto next_ik;
571 ihaka 411 
#endif
1820 ihaka 412 
		sum_r += (xji_r * yjk_r - xji_i * yjk_i);
413 
		sum_i += (xji_r * yjk_i + xji_i * yjk_r);
414 
	    }
415 
	    z[i + k * ncx].r = sum_r;
416 
	    z[i + k * ncx].i = sum_i;
1016 maechler 417 
#ifndef IEEE_754
1820 ihaka 418 
	next_ik:
419 
	    ;
1016 maechler 420 
#endif
1820 ihaka 421 
	}
2 r 422 
}
423 
 
424 
SEXP do_matprod(SEXP call, SEXP op, SEXP args, SEXP rho)
425 
{
1820 ihaka 426 
    int ldx, ldy, nrx, ncx, nry, ncy, mode;
427 
    SEXP x, y, xdims, ydims, ans;
2 r 428 
 
1820 ihaka 429 
    if (!(isNumeric(CAR(args)) || isComplex(CAR(args))) ||
430 
	!(isNumeric(CADR(args)) || isComplex(CADR(args))))
431 
	error("%%*%% requires numeric matrix/vector arguments\n");
2 r 432 
 
1820 ihaka 433 
    x = CAR(args);
434 
    y = CADR(args);
435 
    xdims = getAttrib(x, R_DimSymbol);
436 
    ydims = getAttrib(y, R_DimSymbol);
437 
    ldx = length(xdims);
438 
    ldy = length(ydims);
2 r 439 
 
1820 ihaka 440 
    if (ldx != 2 && ldy != 2) {
1839 ihaka 441 
	if (PRIMVAL(op) == 0) {
1820 ihaka 442 
	    nrx = 1;
443 
	    ncx = LENGTH(x);
2 r 444 
	}
1820 ihaka 445 
	else {
446 
	    nrx = LENGTH(x);
447 
	    ncx = 1;
2 r 448 
	}
1820 ihaka 449 
	nry = LENGTH(y);
450 
	ncy = 1;
451 
    }
452 
    else if (ldx != 2) {
453 
	nry = INTEGER(ydims)[0];
454 
	ncy = INTEGER(ydims)[1];
455 
	nrx = 0;
456 
	ncx = 0;
1839 ihaka 457 
	if (PRIMVAL(op) == 0) {
458 
	    if (LENGTH(x) == nry) {
1820 ihaka 459 
		nrx = 1;
460 
		ncx = LENGTH(x);
2213 maechler 461 
	    }
1820 ihaka 462 
	    else if (LENGTH(x) == ncy) {
463 
		ncx=1;
464 
		nrx=LENGTH(x);
465 
	    }
1839 ihaka 466 
	    if (nry * ncy == 1) {
1820 ihaka 467 
		nrx = LENGTH(x);
468 
		ncx = 1;
469 
	    }
2 r 470 
	}
471 
	else {
1839 ihaka 472 
	    if (LENGTH(x) == nry) {
1820 ihaka 473 
		nrx = LENGTH(x);
474 
		ncx = 1;
475 
	    }
2 r 476 
	}
1820 ihaka 477 
    }
478 
    else if (ldy != 2) {
479 
	nrx = INTEGER(xdims)[0];
480 
	ncx = INTEGER(xdims)[1];
481 
	nry = 0;
482 
	ncy = 0;
1839 ihaka 483 
	if (PRIMVAL(op) == 0) {
1820 ihaka 484 
	    if (LENGTH(y) == ncx) {
485 
		nry = LENGTH(y);
486 
		ncy = 1;
487 
	    }
488 
	    else if (LENGTH(y) == nrx){
489 
		ncy = LENGTH(y);
490 
		nry = 1;
491 
	    }
1839 ihaka 492 
	    if (nrx * ncx == 1) {
1820 ihaka 493 
		ncy = LENGTH(y);
494 
		nry = 1;
495 
	    }
2 r 496 
	}
497 
	else {
1820 ihaka 498 
	    if (LENGTH(y) == nrx) {
499 
		nry = LENGTH(y);
500 
		ncy = 1;
501 
	    }
2 r 502 
	}
1820 ihaka 503 
    }
504 
    else {
505 
	nrx = INTEGER(xdims)[0];
506 
	ncx = INTEGER(xdims)[1];
507 
	nry = INTEGER(ydims)[0];
508 
	ncy = INTEGER(ydims)[1];
509 
    }
2 r 510 
 
1839 ihaka 511 
    if (PRIMVAL(op) == 0) {
512 
	if (ncx != nry)
1820 ihaka 513 
	    errorcall(call, "non-conformable arguments\n");
514 
    }
515 
    else {
1839 ihaka 516 
	if (nrx != nry)
1820 ihaka 517 
	    errorcall(call, "non-conformable arguments\n");
518 
    }
519 
 
1839 ihaka 520 
    if (isComplex(CAR(args)) || isComplex(CADR(args)))
1820 ihaka 521 
	mode = CPLXSXP;
522 
    else
523 
	mode = REALSXP;
524 
    CAR(args) = coerceVector(CAR(args), mode);
525 
    CADR(args) = coerceVector(CADR(args), mode);
526 
 
1839 ihaka 527 
    if (PRIMVAL(op) == 0) {
1820 ihaka 528 
	PROTECT(ans = allocMatrix(mode, nrx, ncy));
1839 ihaka 529 
	if (mode == CPLXSXP)
1820 ihaka 530 
	    cmatprod(COMPLEX(CAR(args)), nrx, ncx,
531 
		     COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));
2 r 532 
	else
1820 ihaka 533 
	    matprod(REAL(CAR(args)), nrx, ncx,
534 
		    REAL(CADR(args)), nry, ncy, REAL(ans));
535 
	PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));
536 
	PROTECT(ydims = getAttrib(CADR(args), R_DimNamesSymbol));
537 
	if (xdims != R_NilValue || ydims != R_NilValue) {
1858 ihaka 538 
	    SEXP dimnames = allocVector(VECSXP, 2);
539 
	    if (xdims != R_NilValue)
540 
		VECTOR(dimnames)[0] = VECTOR(xdims)[0];
541 
	    if (ydims != R_NilValue)
542 
		VECTOR(dimnames)[1] = VECTOR(ydims)[1];
543 
	    setAttrib(ans, R_DimNamesSymbol, dimnames);
2 r 544 
	}
1820 ihaka 545 
    }
546 
    else {
547 
	PROTECT(ans = allocMatrix(mode, ncx, ncy));
1839 ihaka 548 
	if (mode == CPLXSXP)
1820 ihaka 549 
	    ccrossprod(COMPLEX(CAR(args)), nrx, ncx,
550 
		       COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));
551 
	else
552 
	    crossprod(REAL(CAR(args)), nrx, ncx,
553 
		      REAL(CADR(args)), nry, ncy, REAL(ans));
554 
	PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));
555 
	PROTECT(ydims = getAttrib(CADR(args), R_DimNamesSymbol));
556 
	if (xdims != R_NilValue || ydims != R_NilValue) {
1858 ihaka 557 
	    SEXP dimnames = allocVector(VECSXP, 2);
558 
	    if (xdims != R_NilValue)
559 
		VECTOR(dimnames)[0] = VECTOR(xdims)[1];
560 
	    if (ydims != R_NilValue)
561 
		VECTOR(dimnames)[1] = VECTOR(ydims)[1];
562 
	    setAttrib(ans, R_DimNamesSymbol, dimnames);
2 r 563 
	}
1820 ihaka 564 
    }
565 
    UNPROTECT(3);
566 
    return ans;
2 r 567 
}
568 
 
569 
SEXP do_transpose(SEXP call, SEXP op, SEXP args, SEXP rho)
570 
{
1858 ihaka 571 
    SEXP a, r, dims, dimnames, rnames, cnames;
1839 ihaka 572 
    int i, len = 0, ncol=0, nrow=0;
2 r 573 
 
1820 ihaka 574 
    checkArity(op, args);
575 
    a = CAR(args);
2 r 576 
 
1839 ihaka 577 
    if (isVector(a)) {
1820 ihaka 578 
	dims = getAttrib(a, R_DimSymbol);
1858 ihaka 579 
	rnames = R_NilValue;
580 
	cnames = R_NilValue;
1820 ihaka 581 
	switch(length(dims)) {
582 
	case 0:
1858 ihaka 583 
	    nrow = len = length(a);
584 
	    ncol = 1;
585 
	    rnames = getAttrib(a, R_NamesSymbol);
586 
	    break;
1820 ihaka 587 
	case 1:
588 
	    nrow = len = length(a);
589 
	    ncol = 1;
1858 ihaka 590 
	    rnames = getAttrib(a, R_DimNamesSymbol);
591 
	    if (rnames != R_NilValue)
592 
		rnames = VECTOR(rnames)[0];
1820 ihaka 593 
	    break;
594 
	case 2:
595 
	    ncol = ncols(a);
596 
	    nrow = nrows(a);
597 
	    len = length(a);
1858 ihaka 598 
	    dimnames = getAttrib(a, R_DimNamesSymbol);
599 
	    if (dimnames != R_NilValue) {
600 
		rnames = VECTOR(dimnames)[0];
601 
		cnames = VECTOR(dimnames)[1];
602 
	    }
1820 ihaka 603 
	    break;
604 
	default:
605 
	    goto not_matrix;
2 r 606 
	}
1820 ihaka 607 
    }
1858 ihaka 608 
    else
609 
	goto not_matrix;
1820 ihaka 610 
    PROTECT(r = allocVector(TYPEOF(a), len));
611 
    switch (TYPEOF(a)) {
612 
    case LGLSXP:
613 
    case INTSXP:
614 
	for (i = 0; i < len; i++)
615 
	    INTEGER(r)[i] = INTEGER(a)[(i / ncol) + (i % ncol) * nrow];
616 
	break;
617 
    case REALSXP:
618 
	for (i = 0; i < len; i++)
619 
	    REAL(r)[i] = REAL(a)[(i / ncol) + (i % ncol) * nrow];
620 
	break;
621 
    case CPLXSXP:
622 
	for (i = 0; i < len; i++)
623 
	    COMPLEX(r)[i] = COMPLEX(a)[(i / ncol) + (i % ncol) * nrow];
624 
	break;
625 
    case STRSXP:
626 
	for (i = 0; i < len; i++)
627 
	    STRING(r)[i] = STRING(a)[(i / ncol) + (i % ncol) * nrow];
628 
	break;
1858 ihaka 629 
    case VECSXP:
630 
	for (i = 0; i < len; i++)
631 
	    VECTOR(r)[i] = VECTOR(a)[(i / ncol) + (i % ncol) * nrow];
632 
	break;
633 
    default:
634 
	goto not_matrix;
1820 ihaka 635 
    }
1858 ihaka 636 
    PROTECT(dims = allocVector(INTSXP, 2));
1820 ihaka 637 
    INTEGER(dims)[0] = ncol;
638 
    INTEGER(dims)[1] = nrow;
639 
    setAttrib(r, R_DimSymbol, dims);
1858 ihaka 640 
    UNPROTECT(1);
641 
    if(rnames != R_NilValue || cnames != R_NilValue) {
642 
	PROTECT(dimnames = allocVector(VECSXP, 2));
643 
	VECTOR(dimnames)[0] = cnames;
644 
	VECTOR(dimnames)[1] = rnames;
645 
	setAttrib(r, R_DimNamesSymbol, dimnames);
2 r 646 
	UNPROTECT(1);
1820 ihaka 647 
    }
1858 ihaka 648 
    copyMostAttrib(a, r);
1820 ihaka 649 
    UNPROTECT(1);
650 
    return r;
651 
 not_matrix:
652 
    errorcall(call, "argument is not a matrix\n");
653 
    return call;/* never used; just for -Wall */
2 r 654 
}
655 
 
656 
 
1820 ihaka 657 
/* swap works by finding for a index i, the position */
658 
/* in the array with dimensions dims1 in terms of */
659 
/* (i, j, k, l, m, ...); i.e. component-wise position, */
660 
/* then permute these to the order of the array with */
661 
/* dimensions dims2 and work backwards to an integer */
662 
/* offset in this array */
2 r 663 
 
1820 ihaka 664 
static int swap(int ival, SEXP dims1, SEXP dims2, SEXP perm,
665 
		SEXP ind1, SEXP ind2)
2 r 666 
{
1820 ihaka 667 
    int len, t1, i;
668 
    len = length(dims1);
669 
    t1 = ival;
670 
    for (i = 0; i < len; i++) {
671 
	INTEGER(ind1)[i] = t1 % INTEGER(dims1)[i];
672 
	t1 = t1 / INTEGER(dims1)[i];
673 
    }
2 r 674 
 
1820 ihaka 675 
    for (i = 0; i < len; i++)
676 
	INTEGER(ind2)[i] = INTEGER(ind1)[(INTEGER(perm)[i] - 1)];
2 r 677 
 
1820 ihaka 678 
    t1 = INTEGER(ind2)[(len - 1)];
679 
    for (i = (len - 2); i >= 0; i--) {
680 
	t1 *= INTEGER(dims2)[i];
681 
	t1 += INTEGER(ind2)[i];
682 
    }
683 
    return t1;
2 r 684 
}
685 
 
686 
SEXP do_aperm(SEXP call, SEXP op, SEXP args, SEXP rho)
687 
{
1820 ihaka 688 
    SEXP a, perm, r, dimsa, dimsr, ind1, ind2;
689 
    int i, j, len;
2 r 690 
 
1820 ihaka 691 
    checkArity(op, args);
2 r 692 
 
1820 ihaka 693 
    a = CAR(args);
694 
    PROTECT(dimsa = getAttrib(a, R_DimSymbol));
695 
    if (dimsa == R_NilValue)
696 
	error("aperm: invalid first argument, must be an array\n");
2 r 697 
 
1820 ihaka 698 
    PROTECT(perm = coerceVector(CADR(args), INTSXP));
699 
    if (!isVector(perm) || (length(perm) != length(dimsa)))
700 
	error("aperm: invalid second argument, must be a vector\n");
2 r 701 
 
1820 ihaka 702 
    len = length(a);
2 r 703 
 
1820 ihaka 704 
    PROTECT(dimsr = allocVector(INTSXP, length(dimsa)));
705 
    for (i = 0; i < length(dimsa); i++)
706 
	INTEGER(dimsr)[i] = INTEGER(dimsa)[(INTEGER(perm)[i] - 1)];
2 r 707 
 
1820 ihaka 708 
    PROTECT(r = allocVector(TYPEOF(a), len));
709 
    PROTECT(ind1 = allocVector(INTSXP, LENGTH(dimsa)));
710 
    PROTECT(ind2 = allocVector(INTSXP, LENGTH(dimsa)));
2 r 711 
 
1820 ihaka 712 
    switch (TYPEOF(a)) {
713 
    case INTSXP:
714 
    case LGLSXP:
715 
	for (i = 0; i < len; i++) {
716 
	    j = swap(i, dimsa, dimsr, perm, ind1, ind2);
717 
	    INTEGER(r)[j] = INTEGER(a)[i];
2 r 718 
	}
1820 ihaka 719 
	break;
720 
    case REALSXP:
721 
	for (i = 0; i < len; i++) {
722 
	    j = swap(i, dimsa, dimsr, perm, ind1, ind2);
723 
	    REAL(r)[j] = REAL(a)[i];
724 
	}
725 
	break;
726 
    case CPLXSXP:
727 
	for (i = 0; i < len; i++) {
728 
	    j = swap(i, dimsa, dimsr, perm, ind1, ind2);
729 
	    COMPLEX(r)[j] = COMPLEX(a)[i];
730 
	}
731 
	break;
732 
    case STRSXP:
733 
	for (i = 0; i < len; i++) {
734 
	    j = swap(i, dimsa, dimsr, perm, ind1, ind2);
735 
	    STRING(r)[j] = STRING(a)[i];
736 
	}
737 
	break;
1839 ihaka 738 
    case VECSXP:
739 
	for (i = 0; i < len; i++) {
740 
	    j = swap(i, dimsa, dimsr, perm, ind1, ind2);
741 
	    VECTOR(r)[j] = VECTOR(a)[i];
742 
	}
1820 ihaka 743 
    default:
744 
	errorcall(call, "invalid argument\n");
745 
    }
2 r 746 
 
1820 ihaka 747 
    if (INTEGER(CAR(CDDR(args)))[0])
748 
	setAttrib(r, R_DimSymbol, dimsr);
749 
    else
750 
	setAttrib(r, R_DimSymbol, dimsa);
751 
    UNPROTECT(6);
752 
    return r;
2 r 753 
}