The R Project SVN R

Rev

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