The R Project SVN – R – Diff – /branches/djm-source/src/main/array.c


/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
 *  Copyright (C) 1998-2001   The R Development Core Team
 *  Copyright (C) 2002--2004  The R Foundation
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  A copy of the GNU General Public License is available via WWW at
 *  http://www.gnu.org/copyleft/gpl.html.  You can also obtain it by
 *  writing to the Free Software Foundation, Inc., 59 Temple Place,
 *  Suite 330, Boston, MA  02111-1307  USA.
 */
 
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
 
#include <Defn.h>
#include <Rmath.h>
#include <R_ext/RS.h>
#include <R_ext/Applic.h> /* for dgemm */
 
/* "GetRowNames" and "GetColNames" are utility routines which
 * locate and return the row names and column names from the
 * dimnames attribute of a matrix.  They are useful because
 * old versions of R used pair-based lists for dimnames
 * whereas recent versions use vector based lists.
 
 * These are now very old, plus
 * ``When the "dimnames" attribute is
 *   grabbed off an array it is always adjusted to be a vector.''
*/
SEXP GetRowNames(SEXP dimnames)
{
    if (TYPEOF(dimnames) == VECSXP)
	return VECTOR_ELT(dimnames, 0);
    else
	return R_NilValue;
}
 
SEXP GetColNames(SEXP dimnames)
{
    if (TYPEOF(dimnames) == VECSXP)
	return VECTOR_ELT(dimnames, 1);
    else
	return R_NilValue;
}
 
SEXP do_matrix(SEXP call, SEXP op, SEXP args, SEXP rho)
{
    SEXP vals, snr, snc;
    int nr, nc, byrow, lendat;
 
    checkArity(op, args);
    vals = CAR(args);
    snr = CADR(args);
    snc = CADDR(args);
    byrow = asInteger(CADR(CDDR(args)));
 
    /* R wrapper does as.vector
    if (isVector(vals) || isList(vals)) {
	if (length(vals) < 0)  (sic! cannot happen)
	   errorcall(call, "argument has length zero");
    }
    else errorcall(call, "invalid matrix element type"); */
 
    if (!isNumeric(snr) || !isNumeric(snc))
	error("non-numeric matrix extent");
 
    lendat = length(vals);
    nr = asInteger(snr);
    if (nr == NA_INTEGER) /* This is < 0 */
	error("matrix: invalid nrow value (too large or NA)");
    if (nr < 0)
	error("matrix: invalid nrow value (< 0)");
    nc = asInteger(snc);
    if (nc < 0)
	error("matrix: invalid ncol value (< 0)");
    if (nc == NA_INTEGER)
	error("matrix: invalid ncol value (too large or NA)");
    if (nc < 0)
	error("matrix: invalid ncol value (< 0)");
 
    if(lendat > 0 ) {
	if (lendat > 1 && (nr * nc) % lendat != 0) {
	    if (((lendat > nr) && (lendat / nr) * nr != lendat) ||
		((lendat < nr) && (nr / lendat) * lendat != nr))
		warning("data length [%d] is not a sub-multiple or multiple of the number of rows [%d] in matrix", lendat, nr);
	    else if (((lendat > nc) && (lendat / nc) * nc != lendat) ||
		     ((lendat < nc) && (nc / lendat) * lendat != nc))
		warning("data length [%d] is not a sub-multiple or multiple of the number of columns [%d] in matrix", lendat, nc);
	}
	else if ((lendat > 1) && (nr * nc == 0)){
	    warning("data length exceeds size of matrix");
	}
    }
 
    if ((double)nr * (double)nc > INT_MAX)
	error("matrix: too many elements specified");
 
    PROTECT(snr = allocMatrix(TYPEOF(vals), nr, nc));
    if(lendat) {
	if (isVector(vals))
	    copyMatrix(snr, vals, byrow);
	else
	    copyListMatrix(snr, vals, byrow);
    } else if (isVector(vals)) { /* fill with NAs */
	int i, j;
	switch(TYPEOF(vals)) {
	case STRSXP:
	    for (i = 0; i < nr; i++)
		for (j = 0; j < nc; j++)
		    SET_STRING_ELT(snr, i + j * nr, NA_STRING);
	    break;
	case LGLSXP:
	    for (i = 0; i < nr; i++)
		for (j = 0; j < nc; j++)
		    LOGICAL(snr)[i + j * nr] = NA_LOGICAL;
	    break;
	case INTSXP:
	    for (i = 0; i < nr; i++)
		for (j = 0; j < nc; j++)
		    INTEGER(snr)[i + j * nr] = NA_INTEGER;
	    break;
	case REALSXP:
	    for (i = 0; i < nr; i++)
		for (j = 0; j < nc; j++)
		    REAL(snr)[i + j * nr] = NA_REAL;
	    break;
	case CPLXSXP:
	    {
		Rcomplex na_cmplx;
		na_cmplx.r = NA_REAL;
		na_cmplx.i = 0;
		for (i = 0; i < nr; i++)
		    for (j = 0; j < nc; j++)
			COMPLEX(snr)[i + j * nr] = na_cmplx;
	    }
	    break;
	case RAWSXP:
	    for (i = 0; i < nr; i++)
		for (j = 0; j < nc; j++)
		    RAW(snr)[i + j * nr] = 0;
	    break;
	}
    }
    UNPROTECT(1);
    return snr;
}
 
 
SEXP allocMatrix(SEXPTYPE mode, int nrow, int ncol)
{
    SEXP s, t;
    int n;
 
    if (nrow < 0 || ncol < 0)
	error("negative extents to matrix");
    if ((double)nrow * (double)ncol > INT_MAX)
	error("allocMatrix: too many elements specified");
    n = nrow * ncol;
    PROTECT(s = allocVector(mode, n));
    PROTECT(t = allocVector(INTSXP, 2));
    INTEGER(t)[0] = nrow;
    INTEGER(t)[1] = ncol;
    setAttrib(s, R_DimSymbol, t);
    UNPROTECT(2);
    return s;
}
 
 
SEXP allocArray(SEXPTYPE mode, SEXP dims)
{
    SEXP array;
    int i, n;
    double dn;
 
    dn = n = 1;
    for (i = 0; i < LENGTH(dims); i++) {
	dn *= INTEGER(dims)[i];
	if(dn > INT_MAX)
	    error("allocArray: too many elements specified by dims");
	n *= INTEGER(dims)[i];
    }
 
    PROTECT(dims = duplicate(dims));
    PROTECT(array = allocVector(mode, n));
    setAttrib(array, R_DimSymbol, dims);
    UNPROTECT(2);
    return array;
}
 
/* DropDims strips away redundant dimensioning information. */
/* If there is an appropriate dimnames attribute the correct */
/* element is extracted and attached to the vector as a names */
/* attribute.  Note that this function mutates x. */
/* Duplication should occur before this is called. */
 
SEXP DropDims(SEXP x)
{
    SEXP q, dims, dimnames, newnames = R_NilValue;
    int i, n, ndims;
 
    PROTECT(x);
    dims = getAttrib(x, R_DimSymbol);
    dimnames = getAttrib(x, R_DimNamesSymbol);
 
    /* Check that dropping will actually do something. */
    /* (1) Check that there is a "dim" attribute. */
 
    if (dims == R_NilValue) {
	UNPROTECT(1);
	return x;
    }
    ndims = LENGTH(dims);
 
    /* (2) Check whether there are redundant extents */
    n = 0;
    for (i = 0; i < ndims; i++)
	if (INTEGER(dims)[i] != 1) n++;
    if (n == ndims) {
	UNPROTECT(1);
	return x;
    }
 
    if (n <= 1) {
	/* We have reduced to a vector result. */
	if (dimnames != R_NilValue) {
	    n = length(dims);
	    if (TYPEOF(dimnames) == VECSXP) {
		for (i = 0; i < n; i++) {
		    if (INTEGER(dims)[i] != 1) {
			newnames = VECTOR_ELT(dimnames, i);
			break;
		    }
		}
	    }
	    else {
		q = dimnames;
		for (i = 0; i < n; i++) {
		    if (INTEGER(dims)[i] != 1) {
			newnames = CAR(q);
			break;
		    }
		    q = CDR(q);
		}
	    }
	}
	PROTECT(newnames);
	setAttrib(x, R_DimNamesSymbol, R_NilValue);
	setAttrib(x, R_DimSymbol, R_NilValue);
	setAttrib(x, R_NamesSymbol, newnames);
	UNPROTECT(1);
    }
    else {
	/* We have a lower dimensional array. */
	SEXP newdims, dnn, newnamesnames = R_NilValue;
	dnn = getAttrib(dimnames, R_NamesSymbol);
	PROTECT(newdims = allocVector(INTSXP, n));
	for (i = 0, n = 0; i < ndims; i++)
	    if (INTEGER(dims)[i] != 1)
		INTEGER(newdims)[n++] = INTEGER(dims)[i];
	if (!isNull(dimnames)) {
	    int havenames = 0;
	    for (i = 0; i < ndims; i++)
		if (INTEGER(dims)[i] != 1 &&
		    VECTOR_ELT(dimnames, i) != R_NilValue)
		    havenames = 1;
	    if (havenames) {
		PROTECT(newnames = allocVector(VECSXP, n));
		PROTECT(newnamesnames = allocVector(STRSXP, n));
		for (i = 0, n = 0; i < ndims; i++) {
		    if (INTEGER(dims)[i] != 1) {
			if(!isNull(dnn))
			    SET_STRING_ELT(newnamesnames, n,
					   STRING_ELT(dnn, i));
			SET_VECTOR_ELT(newnames, n++, VECTOR_ELT(dimnames, i));
		    }
		}
	    }
	    else dimnames = R_NilValue;
	}
	PROTECT(dimnames);
	setAttrib(x, R_DimNamesSymbol, R_NilValue);
	setAttrib(x, R_DimSymbol, newdims);
	if (dimnames != R_NilValue)
	{
	    if(!isNull(dnn))
		setAttrib(newnames, R_NamesSymbol, newnamesnames);
	    setAttrib(x, R_DimNamesSymbol, newnames);
	    UNPROTECT(2);
	}
	UNPROTECT(2);
    }
    UNPROTECT(1);
    return x;
}
 
SEXP do_drop(SEXP call, SEXP op, SEXP args, SEXP rho)
{
    SEXP x, xdims;
    int i, n, shorten;
 
    checkArity(op, args);
    x = CAR(args);
    if ((xdims = getAttrib(x, R_DimSymbol)) != R_NilValue) {
	n = LENGTH(xdims);
	shorten = 0;
	for (i = 0; i < n; i++)
	    if (INTEGER(xdims)[i] == 1) shorten = 1;
	if (shorten) {
	    if (NAMED(x)) x = duplicate(x);
	    x = DropDims(x);
	}
    }
    return x;
}
 
/* Length of Primitive Objects */
 
SEXP do_length(SEXP call, SEXP op, SEXP args, SEXP rho)
{
    SEXP ans;
    R_len_t len;
 
    if (length(args) != 1)
	error("incorrect number of args to length");
 
    if( isObject(CAR(args)) && DispatchOrEval(call, op, "length", args,
					      rho, &ans, 0, 1))
      return(ans);
 
    ans = allocVector(INTSXP, 1);
    len = length(CAR(args));
    INTEGER(ans)[0] = (len < INT_MAX) ? len : NA_INTEGER;
    return ans;
}
 
 
SEXP do_rowscols(SEXP call, SEXP op, SEXP args, SEXP rho)
{
    SEXP ans;
    int i, j, nr, nc;
 
    if (length(args) != 1)
	error("incorrect number of args to row/col");
    if (!isMatrix(CAR(args)))
	error("a matrix is required as arg to row/col");
 
    nr = nrows(CAR(args));
    nc = ncols(CAR(args));
 
    ans = allocMatrix(INTSXP, nr, nc);
 
    switch (PRIMVAL(op)) {
    case 1:
	for (i = 0; i < nr; i++)
	    for (j = 0; j < nc; j++)
		INTEGER(ans)[i + j * nr] = i + 1;
	break;
    case 2:
	for (i = 0; i < nr; i++)
	    for (j = 0; j < nc; j++)
		INTEGER(ans)[i + j * nr] = j + 1;
	break;
    }
    return ans;
}
 
static void matprod(double *x, int nrx, int ncx,
		    double *y, int nry, int ncy, double *z)
{
    char *transa = "N", *transb = "N";
    int i,  j, k;
    double one = 1.0, zero = 0.0, sum;
    Rboolean have_na = FALSE;
 
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
	/* Don't trust the BLAS to handle NA/NaNs correctly: PR#4582
	 * The test is only O(n) here
	 */
	for (i = 0; i < nrx*ncx; i++)
	    if (ISNAN(x[i])) {have_na = TRUE; break;}
	if (!have_na) 
	    for (i = 0; i < nry*ncy; i++)
		if (ISNAN(y[i])) {have_na = TRUE; break;}
	if (have_na) {
	    for (i = 0; i < nrx; i++)
		for (k = 0; k < ncy; k++) {
		    sum = 0.0;
		    for (j = 0; j < ncx; j++)
			sum += x[i + j * nrx] * y[j + k * nry];
		    z[i + k * nrx] = sum;
		}
	} else
	    F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
			    x, &nrx, y, &nry, &zero, z, &nrx);
    } else /* zero-extent operations should return zeroes */
	for(i = 0; i < nrx*ncy; i++) z[i] = 0;
}
 
#ifdef HAVE_DOUBLE_COMPLEX
/* ZGEMM - perform one of the matrix-matrix operations    */
/* C := alpha*op( A )*op( B ) + beta*C */
extern void
F77_NAME(zgemm)(const char *transa, const char *transb, const int *m,
		const int *n, const int *k, const Rcomplex *alpha,
		const Rcomplex *a, const int *lda,
		const Rcomplex *b, const int *ldb,
		const Rcomplex *beta, Rcomplex *c, const int *ldc);
#endif
 
static void cmatprod(Rcomplex *x, int nrx, int ncx,
		     Rcomplex *y, int nry, int ncy, Rcomplex *z)
{
#ifdef HAVE_DOUBLE_COMPLEX
    char *transa = "N", *transb = "N";
    int i;
    Rcomplex one, zero;
 
    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
        F77_CALL(zgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
			x, &nrx, y, &nry, &zero, z, &nrx);
    } else { /* zero-extent operations should return zeroes */
	for(i = 0; i < nrx*ncy; i++) z[i].r = z[i].i = 0;
    }
#else
    int i, j, k;
    double xij_r, xij_i, yjk_r, yjk_i, sum_i, sum_r;
 
    for (i = 0; i < nrx; i++)
	for (k = 0; k < ncy; k++) {
	    z[i + k * nrx].r = NA_REAL;
	    z[i + k * nrx].i = NA_REAL;
	    sum_r = 0.0;
	    sum_i = 0.0;
	    for (j = 0; j < ncx; j++) {
		xij_r = x[i + j * nrx].r;
		xij_i = x[i + j * nrx].i;
		yjk_r = y[j + k * nry].r;
		yjk_i = y[j + k * nry].i;
		if (ISNAN(xij_r) || ISNAN(xij_i)
		    || ISNAN(yjk_r) || ISNAN(yjk_i))
		    goto next_ik;
		sum_r += (xij_r * yjk_r - xij_i * yjk_i);
		sum_i += (xij_r * yjk_i + xij_i * yjk_r);
	    }
	    z[i + k * nrx].r = sum_r;
	    z[i + k * nrx].i = sum_i;
	next_ik:
	    ;
	}
#endif
}
 
static void symcrossprod(double *x, int nr, int nc, double *z)
{
    char *trans = "T", *uplo = "U";
    double one = 1.0, zero = 0.0;
    int i, j;
    if (nr > 0 && nc > 0) {
        F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc);
	for (i = 1; i < nc; i++)
	    for (j = 0; j < i; j++) z[i + nc *j] = z[j + nc * i];
    } else { /* zero-extent operations should return zeroes */
	for(i = 0; i < nc*nc; i++) z[i] = 0;
    }
 
}
 
static void crossprod(double *x, int nrx, int ncx,
		      double *y, int nry, int ncy, double *z)
{
    char *transa = "T", *transb = "N";
    double one = 1.0, zero = 0.0;
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
        F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
			x, &nrx, y, &nry, &zero, z, &ncx);
    }
    else { /* zero-extent operations should return zeroes */
	int i;
	for(i = 0; i < ncx*ncy; i++) z[i] = 0;
    }
}
 
static void ccrossprod(Rcomplex *x, int nrx, int ncx,
		       Rcomplex *y, int nry, int ncy, Rcomplex *z)
{
    char *transa = "T", *transb = "N";
    Rcomplex one, zero;
 
    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
        F77_CALL(zgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
			x, &nrx, y, &nry, &zero, z, &ncx);
    }
    else { /* zero-extent operations should return zeroes */
	int i;
	for(i = 0; i < ncx*ncy; i++) z[i].r = z[i].i = 0;
    }
}
/* "%*%" (op = 0)  or  crossprod (op = 1) : */
SEXP do_matprod(SEXP call, SEXP op, SEXP args, SEXP rho)
{
    int ldx, ldy, nrx, ncx, nry, ncy, mode;
    SEXP x = CAR(args), y = CADR(args), xdims, ydims, ans;
    Rboolean sym;
 
    if(R_has_methods(op)) {
      SEXP value;
      value = R_possible_dispatch(call, op, args, rho);
      if(value) return value;
    }
 
    sym = isNull(y);
    if (sym && (PRIMVAL(op) == 1)) y = x;
    if ( !(isNumeric(x) || isComplex(x)) || !(isNumeric(y) || isComplex(y)) )
	errorcall(call, "requires numeric matrix/vector arguments");
 
    xdims = getAttrib(x, R_DimSymbol);
    ydims = getAttrib(y, R_DimSymbol);
    ldx = length(xdims);
    ldy = length(ydims);
 
    if (ldx != 2 && ldy != 2) {		/* x and y non-matrices */
	if (PRIMVAL(op) == 0) {
	    nrx = 1;
	    ncx = LENGTH(x);
	}
	else {
	    nrx = LENGTH(x);
	    ncx = 1;
	}
	nry = LENGTH(y);
	ncy = 1;
    }
    else if (ldx != 2) {		/* x not a matrix */
	nry = INTEGER(ydims)[0];
	ncy = INTEGER(ydims)[1];
	nrx = 0;
	ncx = 0;
	if (PRIMVAL(op) == 0) {
	    if (LENGTH(x) == nry) {	/* x as row vector */
		nrx = 1;
		ncx = LENGTH(x);
	    }
	    else if (nry == 1) {	/* x as col vector */
		nrx = LENGTH(x);
		ncx = 1;
	    }
	}
	else {
	    if (LENGTH(x) == nry) {	/* x is a row vector */
		nrx = LENGTH(x);
		ncx = 1;
	    }
	}
    }
    else if (ldy != 2) {		/* y not a matrix */
	nrx = INTEGER(xdims)[0];
	ncx = INTEGER(xdims)[1];
	nry = 0;
	ncy = 0;
	if (PRIMVAL(op) == 0) {
	    if (LENGTH(y) == ncx) {	/* y as col vector */
		nry = LENGTH(y);
		ncy = 1;
	    }
	    else if (ncx == 1) {	/* y as row vector */
		nry = 1;
		ncy = LENGTH(y);
	    }
	}
	else {
	    if (LENGTH(y) == nrx) {	/* y is a row vector */
		nry = LENGTH(y);
		ncy = 1;
	    }
	}
    }
    else {				/* x and y matrices */
	nrx = INTEGER(xdims)[0];
	ncx = INTEGER(xdims)[1];
	nry = INTEGER(ydims)[0];
	ncy = INTEGER(ydims)[1];
    }
 
    if (PRIMVAL(op) == 0) {
	if (ncx != nry)
	    errorcall(call, "non-conformable arguments");
    }
    else {
	if (nrx != nry)
	    errorcall(call, "non-conformable arguments");
    }
 
    if (isComplex(CAR(args)) || isComplex(CADR(args)))
	mode = CPLXSXP;
    else
	mode = REALSXP;
    SETCAR(args, coerceVector(CAR(args), mode));
    SETCADR(args, coerceVector(CADR(args), mode));
 
    if (PRIMVAL(op) == 0) {		       	/* op == 0 : matprod() */
 
	PROTECT(ans = allocMatrix(mode, nrx, ncy));
	if (mode == CPLXSXP)
	    cmatprod(COMPLEX(CAR(args)), nrx, ncx,
		     COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));
	else
	    matprod(REAL(CAR(args)), nrx, ncx,
		    REAL(CADR(args)), nry, ncy, REAL(ans));
 
	PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));
	PROTECT(ydims = getAttrib(CADR(args), R_DimNamesSymbol));
 
	if (xdims != R_NilValue || ydims != R_NilValue) {
	    SEXP dimnames, dimnamesnames, dn;
	    PROTECT(dimnames = allocVector(VECSXP, 2));
	    PROTECT(dimnamesnames = allocVector(STRSXP, 2));
	    if (xdims != R_NilValue) {
		if (ldx == 2 || ncx ==1) {
		    dn = getAttrib(xdims, R_NamesSymbol);
		    SET_VECTOR_ELT(dimnames, 0, VECTOR_ELT(xdims, 0));
		    if(!isNull(dn))
			SET_STRING_ELT(dimnamesnames, 0, STRING_ELT(dn, 0));
		}
	    }
	    if (ydims != R_NilValue) {
		if (ldy == 2 ){
		    dn = getAttrib(ydims, R_NamesSymbol);
		    SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 1));
		    if(!isNull(dn))
			SET_STRING_ELT(dimnamesnames, 1, STRING_ELT(dn, 1));
		} else if (nry == 1) {
		    dn = getAttrib(ydims, R_NamesSymbol);
		    SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 0));
		    if(!isNull(dn))
			SET_STRING_ELT(dimnamesnames, 1, STRING_ELT(dn, 0));
		}
	    }
	    setAttrib(dimnames, R_NamesSymbol, dimnamesnames);
	    setAttrib(ans, R_DimNamesSymbol, dimnames);
	    UNPROTECT(2);
	}
    }
 
    else {					/* op == 1: crossprod() */
 
	PROTECT(ans = allocMatrix(mode, ncx, ncy));
	if (mode == CPLXSXP)
	    if(sym)
		ccrossprod(COMPLEX(CAR(args)), nrx, ncx,
			   COMPLEX(CAR(args)), nry, ncy, COMPLEX(ans));
	    else
		ccrossprod(COMPLEX(CAR(args)), nrx, ncx,
			   COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans));
	else {
	    if(sym)
		symcrossprod(REAL(CAR(args)), nrx, ncx, REAL(ans));
	    else
		crossprod(REAL(CAR(args)), nrx, ncx,
			  REAL(CADR(args)), nry, ncy, REAL(ans));
	}
 
	PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));
	if (sym)
	    PROTECT(ydims = xdims);
	else
	    PROTECT(ydims = getAttrib(CADR(args), R_DimNamesSymbol));
 
	if (xdims != R_NilValue || ydims != R_NilValue) {
	    SEXP dimnames, dimnamesnames, dnx=R_NilValue, dny=R_NilValue;
 
	    /* allocate dimnames and dimnamesnames */
 
	    PROTECT(dimnames = allocVector(VECSXP, 2));
	    PROTECT(dimnamesnames = allocVector(STRSXP, 2));
 
            /* There was a bug here.  The second element of a */
            /* dimnames list was being accessed for a 1-d array. */
            /* I have just excluded the use of dimnames in this */
            /* case. - ihaka Sep 30, 2003. */
 
	    if (xdims != R_NilValue) {
                if (ldx == 2) {
		    dnx = getAttrib(xdims, R_NamesSymbol);
		    SET_VECTOR_ELT(dimnames, 0, VECTOR_ELT(xdims, 1));
		    if(!isNull(dnx))
			SET_STRING_ELT(dimnamesnames, 0, STRING_ELT(dnx, 1));
		}
	    }
 
	    if (ydims != R_NilValue) {
                if (ldy == 2) {
		    dny = getAttrib(ydims, R_NamesSymbol);
		    SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 1));
		    if(!isNull(dny))
			SET_STRING_ELT(dimnamesnames, 1, STRING_ELT(dny, 1));
		}
	    }
 
            /* We sometimes attach a dimnames attribute */
            /* whose elements are all NULL ... */
            /* Thus is ugly but causes no real damage. */
	    
	    if (!isNull(dnx) || !isNull(dny))
		setAttrib(dimnames, R_NamesSymbol, dimnamesnames);
	    setAttrib(ans, R_DimNamesSymbol, dimnames);
	    UNPROTECT(2);
	}
    }
    UNPROTECT(3);
    return ans;
}
 
SEXP do_transpose(SEXP call, SEXP op, SEXP args, SEXP rho)
{
    SEXP a, r, dims, dimnames, dimnamesnames=R_NilValue,
	ndimnamesnames, rnames, cnames;
    int i, len = 0, ncol=0, nrow=0;
 
    checkArity(op, args);
    a = CAR(args);
 
    if (isVector(a)) {
	dims = getAttrib(a, R_DimSymbol);
	rnames = R_NilValue;
	cnames = R_NilValue;
	switch(length(dims)) {
	case 0:
	    nrow = len = length(a);
	    ncol = 1;
	    rnames = getAttrib(a, R_NamesSymbol);
	    break;
	case 1:
	    nrow = len = length(a);
	    ncol = 1;
	    rnames = getAttrib(a, R_DimNamesSymbol);
	    if (rnames != R_NilValue)
		rnames = VECTOR_ELT(rnames, 0);
	    break;
	case 2:
	    ncol = ncols(a);
	    nrow = nrows(a);
	    len = length(a);
	    dimnames = getAttrib(a, R_DimNamesSymbol);
	    if (dimnames != R_NilValue) {
		rnames = VECTOR_ELT(dimnames, 0);
		cnames = VECTOR_ELT(dimnames, 1);
		dimnamesnames = getAttrib(dimnames, R_NamesSymbol);
	    }
	    break;
	default:
	    goto not_matrix;
	}
    }
    else
	goto not_matrix;
    PROTECT(r = allocVector(TYPEOF(a), len));
    switch (TYPEOF(a)) {
    case LGLSXP:
    case INTSXP:
	for (i = 0; i < len; i++)
	    INTEGER(r)[i] = INTEGER(a)[(i / ncol) + (i % ncol) * nrow];
	break;
    case REALSXP:
	for (i = 0; i < len; i++)
	    REAL(r)[i] = REAL(a)[(i / ncol) + (i % ncol) * nrow];
	break;
    case CPLXSXP:
	for (i = 0; i < len; i++)
	    COMPLEX(r)[i] = COMPLEX(a)[(i / ncol) + (i % ncol) * nrow];
	break;
    case STRSXP:
	for (i = 0; i < len; i++)
	    SET_STRING_ELT(r, i,
			   STRING_ELT(a, (i / ncol) + (i % ncol) * nrow));
	break;
    case VECSXP:
	for (i = 0; i < len; i++)
	    SET_VECTOR_ELT(r, i,
			   VECTOR_ELT(a, (i / ncol) + (i % ncol) * nrow));
	break;
    case RAWSXP:
	for (i = 0; i < len; i++)
	    RAW(r)[i] = RAW(a)[(i / ncol) + (i % ncol) * nrow];
	break;
    default:
	goto not_matrix;
    }
    PROTECT(dims = allocVector(INTSXP, 2));
    INTEGER(dims)[0] = ncol;
    INTEGER(dims)[1] = nrow;
    setAttrib(r, R_DimSymbol, dims);
    UNPROTECT(1);
    if(rnames != R_NilValue || cnames != R_NilValue) {
	PROTECT(dimnames = allocVector(VECSXP, 2));
	SET_VECTOR_ELT(dimnames, 0, cnames);
	SET_VECTOR_ELT(dimnames, 1, rnames);
	if(!isNull(dimnamesnames)) {
	    PROTECT(ndimnamesnames = allocVector(VECSXP, 2));
	    SET_STRING_ELT(ndimnamesnames, 1, STRING_ELT(dimnamesnames, 0));
	    SET_STRING_ELT(ndimnamesnames, 0, STRING_ELT(dimnamesnames, 1));
	    setAttrib(dimnames, R_NamesSymbol, ndimnamesnames);
	    UNPROTECT(1);
	}
	setAttrib(r, R_DimNamesSymbol, dimnames);
	UNPROTECT(1);
    }
    copyMostAttrib(a, r);
    UNPROTECT(1);
    return r;
 not_matrix:
    errorcall(call, "argument is not a matrix");
    return call;/* never used; just for -Wall */
}
 
/*
 New version of aperm, using strides for speed.
 Jonathan Rougier <J.C.Rougier@durham.ac.uk>
 
 v1.0 30.01.01
 
 M.Maechler : expanded	all ../include/Rdefines.h macros
 */
 
/* this increments iip and sets j using strides */
 
#define CLICKJ						\
    for (itmp=0; itmp<n; itmp++)			\
	if (iip[itmp] == INTEGER(dimsr)[itmp]-1)	\
	    iip[itmp] = 0;				\
	else {						\
	    iip[itmp]++;				\
	    break;					\
	}						\
    for (j=0, itmp=0; itmp<n; itmp++)			\
	j += iip[itmp] * stride[itmp];
 
/* aperm (a, perm, resize = TRUE) */
SEXP do_aperm(SEXP call, SEXP op, SEXP args, SEXP rho)
{
    SEXP a, perm, resize, r, dimsa, dimsr, dna;
    int i, j, n, len, itmp;
    int *pp, *iip, *stride;
    char *vmax;
 
    checkArity(op, args);
 
    a = CAR(args);
    if (!isArray(a))
	errorcall(call,"invalid first argument, must be an array");
 
    PROTECT(dimsa = getAttrib(a, R_DimSymbol));
    n = LENGTH(dimsa);
 
    /* check the permutation */
 
    PROTECT(perm = coerceVector(CADR(args), INTSXP));
    vmax = vmaxget();
    pp = (int *) R_alloc(n, sizeof(int));
    if (length(perm) == 0) {
	for (i=0; i<n; i++)
	    pp[i] = n-1-i;
    } else if (length(perm) == n) {
	for (i=0; i<n; i++)
	    pp[i] = INTEGER(perm)[i] - 1; /* no offset! */
    } else
	errorcall(call, "`perm' is of wrong length");
 
    iip = (int *) R_alloc(n, sizeof(int));
    for (i=0; i<n; iip[i++] = 0);
    for (i=0; i<n; i++)
	if (pp[i] >= 0 && pp[i] < n)
	    iip[pp[i]]++;
	else
	    errorcall(call, "value of out range in `perm'");
    for (i=0; i<n; i++)
	if (iip[i]==0)
	    errorcall(call, "invalid permutation (`perm')");
 
    /* create the stride object and permute */
 
    stride = (int *) R_alloc(n, sizeof(int));
 
    for (iip[0] = 1, i = 1; i<n; i++)
	iip[i] = iip[i-1] * INTEGER(dimsa)[i-1];
 
    for (i=0; i<n; i++)
	stride[i] = iip[pp[i]];
 
    /* also need to have the dimensions of r */
 
    PROTECT(dimsr = allocVector(INTSXP,n));
    for (i=0; i<n; i++)
	INTEGER(dimsr)[i] = INTEGER(dimsa)[pp[i]];
 
    /* and away we go! iip will hold the incrementer */
 
    len = LENGTH(a);
    len = length(a);
    PROTECT(r = allocVector(TYPEOF(a), len));
 
    for (i=0; i<n; iip[i++] = 0);
 
    switch (TYPEOF(a)) {
 
    case INTSXP:
    case LGLSXP:
	for (j=0, i=0; i<len; i++) {
	    INTEGER(r)[i] = INTEGER(a)[j];
	    CLICKJ;
	}
	break;
 
    case REALSXP:
	for (j=0, i=0; i<len; i++) {
	    REAL(r)[i] = REAL(a)[j];
	    CLICKJ;
	}
	break;
 
    case CPLXSXP:
	for (j=0, i=0; i<len; i++) {
	    COMPLEX(r)[i].r = COMPLEX(a)[j].r;
	    COMPLEX(r)[i].i = COMPLEX(a)[j].i;
	    CLICKJ;
	}
	break;
 
    case STRSXP:
	for (j=0, i=0; i<len; i++) {
	    SET_STRING_ELT(r, i, STRING_ELT(a, j));
	    CLICKJ;
	}
	break;
 
    case VECSXP:
	for (j=0, i=0; i<len; i++) {
	    SET_VECTOR_ELT(r, i, VECTOR_ELT(a, j));
	    CLICKJ;
	}
	break;
 
    case RAWSXP:
	for (j=0, i=0; i<len; i++) {
	    RAW(r)[i] = RAW(a)[j];
	    CLICKJ;
	}
	break;
 
    default:
	errorcall(call, "unsupported type of array");
    }
 
    /* handle the resize */
    PROTECT(resize = coerceVector(CADDR(args), INTSXP));
    if (LOGICAL(resize)[0])
	setAttrib(r, R_DimSymbol, dimsr);
    else
	setAttrib(r, R_DimSymbol, dimsa);
 
    /* and handle the dimnames */
 
    PROTECT(dna = getAttrib(a, R_DimNamesSymbol));
 
    if (LOGICAL(resize)[0] && dna != R_NilValue) {
 
	SEXP dnna, dnr, dnnr;
 
	PROTECT(dnna = getAttrib(dna, R_NamesSymbol));
	PROTECT(dnnr = allocVector(STRSXP,n));
	PROTECT(dnr  = allocVector(VECSXP,n));
 
	for (i=0; i<n; i++) {
	    SET_VECTOR_ELT(dnr, i, VECTOR_ELT(dna, pp[i]));
	    if (dnna != R_NilValue)
		SET_STRING_ELT(dnnr, i, STRING_ELT(dnna, pp[i]));
	}
 
	if (dnna != R_NilValue)
	    setAttrib(dnr, R_NamesSymbol, dnnr);
	setAttrib(r, R_DimNamesSymbol, dnr);
	UNPROTECT(3); /* dnna, dnr, dnnr */
    }
    /* free temporary memory */
    vmaxset(vmax);
 
    UNPROTECT(6); /* dimsa, perm, r, dimsr, resize, dna */
    return r;
}
 
/* colSums(x, n, p, na.rm) and friends */
SEXP do_colsum(SEXP call, SEXP op, SEXP args, SEXP rho)
{
    SEXP x, ans = R_NilValue;
    int OP, n, p, cnt = 0, i, j, type;
    Rboolean NaRm, keepNA;
    int *ix;
    double *rx, sum = 0.0;
 
    checkArity(op, args);
    x = CAR(args); args = CDR(args);
    n = asInteger(CAR(args)); args = CDR(args);
    p = asInteger(CAR(args)); args = CDR(args);
    NaRm = asLogical(CAR(args));
    if (n == NA_INTEGER || n <= 0)
	errorcall(call, "invalid value of n");
    if (p == NA_INTEGER || p <= 0)
	errorcall(call, "invalid value of p");
    if (NaRm == NA_LOGICAL) errorcall(call, "invalid value of na.rm");
    keepNA = !NaRm;
 
    OP = PRIMVAL(op);
    switch (type = TYPEOF(x)) {
    case LGLSXP: break;
    case INTSXP: break;
    case REALSXP: break;
    default:
	errorcall(call, "`x' must be numeric");
    }
 
    if (OP == 0 || OP == 1) { /* columns */
	cnt = n;
	PROTECT(ans = allocVector(REALSXP, p));
	for (j = 0; j < p; j++) {
	    switch (type) {
	    case REALSXP:
		rx = REAL(x) + n*j;
		if (keepNA)
		    for (sum = 0., i = 0; i < n; i++) sum += *rx++;
		else {
		    for (cnt = 0, sum = 0., i = 0; i < n; i++, rx++)
			if (!ISNAN(*rx)) {cnt++; sum += *rx;}
			else if (keepNA) {sum = NA_REAL; break;}
		}
		break;
	    case INTSXP:
		ix = INTEGER(x) + n*j;
		for (cnt = 0, sum = 0., i = 0; i < n; i++, ix++)
		    if (*ix != NA_INTEGER) {cnt++; sum += *ix;}
		    else if (keepNA) {sum = NA_REAL; break;}
		break;
	    case LGLSXP:
		ix = LOGICAL(x) + n*j;
		for (cnt = 0, sum = 0., i = 0; i < n; i++, ix++)
		    if (*ix != NA_LOGICAL) {cnt++; sum += *ix;}
		    else if (keepNA) {sum = NA_REAL; break;}
		break;
	    }
	    if (OP == 1) {
		if (cnt > 0) sum /= cnt; else sum = NA_REAL;
	    }
	    REAL(ans)[j] = sum;
	}
    }
 
    if (OP == 2 || OP == 3) { /* rows */
	cnt = p;
	PROTECT(ans = allocVector(REALSXP, n));
 
	/* reverse summation order to improve cache hits */
	if (type == REALSXP) {
	    double *rans = REAL(ans), *ra = rans, *cnt = NULL, *c;
	    rx = REAL(x);
	    if (!keepNA && OP == 3) cnt = Calloc(n, double);
	    for (ra = rans, i = 0; i < n; i++) *ra++ = 0.0;
	    for (j = 0; j < p; j++) {
		ra = rans;
		if (keepNA)
		    for (i = 0; i < n; i++) *ra++ += *rx++;
		else
		    for (c = cnt, i = 0; i < n; i++, ra++, rx++, c++)
			if (!ISNAN(*rx)) {
			    *ra += *rx;
			    if (OP == 3) (*c)++;
			}
	    }
	    if (OP == 3) {
		if (keepNA)
		    for (ra = rans, i = 0; i < n; i++)
			*ra++ /= p;
		else {
		    for (ra = rans, c = cnt, i = 0; i < n; i++, c++)
			if (*c > 0) *ra++ /= *c; else *ra++ = NA_REAL;
		    Free(cnt);
		}
	    }
	    UNPROTECT(1);
	    return ans;
	}
 
	for (i = 0; i < n; i++) {
	    switch (type) {
	    case REALSXP: /* this cannot be reached */
		rx = REAL(x) + i;
		if (keepNA)
		    for (sum = 0., j = 0; j < p; j++, rx += n) sum += *rx;
		else {
		    for (cnt = 0, sum = 0., j = 0; j < p; j++, rx += n)
			if (!ISNAN(*rx)) {cnt++; sum += *rx;}
			else if (keepNA) {sum = NA_REAL; break;}
		}
		break;
	    case INTSXP:
		ix = INTEGER(x) + i;
		for (cnt = 0, sum = 0., j = 0; j < p; j++, ix += n)
		    if (*ix != NA_INTEGER) {cnt++; sum += *ix;}
		    else if (keepNA) {sum = NA_REAL; break;}
		break;
	    case LGLSXP:
		ix = LOGICAL(x) + i;
		for (cnt = 0, sum = 0., j = 0; j < p; j++, ix += n)
		    if (*ix != NA_LOGICAL) {cnt++; sum += *ix;}
		    else if (keepNA) {sum = NA_REAL; break;}
		break;
	    }
	    if (OP == 3) {
		if (cnt > 0) sum /= cnt; else sum = NA_REAL;
	    }
	    REAL(ans)[i] = sum;
	}
    }
 
    UNPROTECT(1);
    return ans;
}

Rev 29864	Rev 30691
1	`/*`	1	`/*`
2	`* R : A Computer Language for Statistical Data Analysis`	2	`* R : A Computer Language for Statistical Data Analysis`
3	`* Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka`	3	`* Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka`
4	`* Copyright (C) 1998-2001 The R Development Core Team`	4	`* Copyright (C) 1998-2001 The R Development Core Team`
5	`* Copyright (C) 2002--2004 The R Foundation`	5	`* Copyright (C) 2002--2004 The R Foundation`
6	`*`	6	`*`
7	`* This program is free software; you can redistribute it and/or modify`	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`	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`	9	`* the Free Software Foundation; either version 2 of the License, or`
10	`* (at your option) any later version.`	10	`* (at your option) any later version.`
11	`*`	11	`*`
12	`* This program is distributed in the hope that it will be useful,`	12	`* This program is distributed in the hope that it will be useful,`
13	`* but WITHOUT ANY WARRANTY; without even the implied warranty of`	13	`* but WITHOUT ANY WARRANTY; without even the implied warranty of`
14	`* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the`	14	`* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the`
15	`* GNU General Public License for more details.`	15	`* GNU General Public License for more details.`
16	`*`	16	`*`
17	`* A copy of the GNU General Public License is available via WWW at`	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`	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,`	19	`* writing to the Free Software Foundation, Inc., 59 Temple Place,`
20	`* Suite 330, Boston, MA 02111-1307 USA.`	20	`* Suite 330, Boston, MA 02111-1307 USA.`
21	`*/`	21	`*/`
22		22
23	`#ifdef HAVE_CONFIG_H`	23	`#ifdef HAVE_CONFIG_H`
24	`#include <config.h>`	24	`#include <config.h>`
25	`#endif`	25	`#endif`
26		26
27	`#include <Defn.h>`	27	`#include <Defn.h>`
28	`#include <Rmath.h>`	28	`#include <Rmath.h>`
29	`#include <R_ext/RS.h>`	29	`#include <R_ext/RS.h>`
30	`#include <R_ext/Applic.h> /* for dgemm */`	30	`#include <R_ext/Applic.h> /* for dgemm */`
31		31
32	`/* "GetRowNames" and "GetColNames" are utility routines which`	32	`/* "GetRowNames" and "GetColNames" are utility routines which`
33	`* locate and return the row names and column names from the`	33	`* locate and return the row names and column names from the`
34	`* dimnames attribute of a matrix. They are useful because`	34	`* dimnames attribute of a matrix. They are useful because`
35	`* old versions of R used pair-based lists for dimnames`	35	`* old versions of R used pair-based lists for dimnames`
36	`* whereas recent versions use vector based lists.`	36	`* whereas recent versions use vector based lists.`
37		37
38	`* These are now very old, plus`	38	`* These are now very old, plus`
39	* ``When the "dimnames" attribute is	39	* ``When the "dimnames" attribute is
40	`* grabbed off an array it is always adjusted to be a vector.''`	40	`* grabbed off an array it is always adjusted to be a vector.''`
41	`*/`	41	`*/`
42	`SEXP GetRowNames(SEXP dimnames)`	42	`SEXP GetRowNames(SEXP dimnames)`
43	`{`	43	`{`
44	`if (TYPEOF(dimnames) == VECSXP)`	44	`if (TYPEOF(dimnames) == VECSXP)`
45	`return VECTOR_ELT(dimnames, 0);`	45	`return VECTOR_ELT(dimnames, 0);`
46	`else`	46	`else`
47	`return R_NilValue;`	47	`return R_NilValue;`
48	`}`	48	`}`
49		49
50	`SEXP GetColNames(SEXP dimnames)`	50	`SEXP GetColNames(SEXP dimnames)`
51	`{`	51	`{`
52	`if (TYPEOF(dimnames) == VECSXP)`	52	`if (TYPEOF(dimnames) == VECSXP)`
53	`return VECTOR_ELT(dimnames, 1);`	53	`return VECTOR_ELT(dimnames, 1);`
54	`else`	54	`else`
55	`return R_NilValue;`	55	`return R_NilValue;`
56	`}`	56	`}`
57		57
58	`SEXP do_matrix(SEXP call, SEXP op, SEXP args, SEXP rho)`	58	`SEXP do_matrix(SEXP call, SEXP op, SEXP args, SEXP rho)`
59	`{`	59	`{`
60	`SEXP vals, snr, snc;`	60	`SEXP vals, snr, snc;`
61	`int nr, nc, byrow, lendat;`	61	`int nr, nc, byrow, lendat;`
62		62
63	`checkArity(op, args);`	63	`checkArity(op, args);`
64	`vals = CAR(args);`	64	`vals = CAR(args);`
65	`snr = CADR(args);`	65	`snr = CADR(args);`
66	`snc = CADDR(args);`	66	`snc = CADDR(args);`
67	`byrow = asInteger(CADR(CDDR(args)));`	67	`byrow = asInteger(CADR(CDDR(args)));`
68		68
69	`/* R wrapper does as.vector`	69	`/* R wrapper does as.vector`
70	`if (isVector(vals) \|\| isList(vals)) {`	70	`if (isVector(vals) \|\| isList(vals)) {`
71	`if (length(vals) < 0) (sic! cannot happen)`	71	`if (length(vals) < 0) (sic! cannot happen)`
72	`errorcall(call, "argument has length zero");`	72	`errorcall(call, "argument has length zero");`
73	`}`	73	`}`
74	`else errorcall(call, "invalid matrix element type"); */`	74	`else errorcall(call, "invalid matrix element type"); */`
75		75
76	`if (!isNumeric(snr) \|\| !isNumeric(snc))`	76	`if (!isNumeric(snr) \|\| !isNumeric(snc))`
77	`error("non-numeric matrix extent");`	77	`error("non-numeric matrix extent");`
78		78
79	`lendat = length(vals);`	79	`lendat = length(vals);`
80	`nr = asInteger(snr);`	80	`nr = asInteger(snr);`
81	`if (nr == NA_INTEGER) /* This is < 0 */`	81	`if (nr == NA_INTEGER) /* This is < 0 */`
82	`error("matrix: invalid nrow value (too large or NA)");`	82	`error("matrix: invalid nrow value (too large or NA)");`
83	`if (nr < 0)`	83	`if (nr < 0)`
84	`error("matrix: invalid nrow value (< 0)");`	84	`error("matrix: invalid nrow value (< 0)");`
85	`nc = asInteger(snc);`	85	`nc = asInteger(snc);`
86	`if (nc < 0)`	86	`if (nc < 0)`
87	`error("matrix: invalid ncol value (< 0)");`	87	`error("matrix: invalid ncol value (< 0)");`
88	`if (nc == NA_INTEGER)`	88	`if (nc == NA_INTEGER)`
89	`error("matrix: invalid ncol value (too large or NA)");`	89	`error("matrix: invalid ncol value (too large or NA)");`
90	`if (nc < 0)`	90	`if (nc < 0)`
91	`error("matrix: invalid ncol value (< 0)");`	91	`error("matrix: invalid ncol value (< 0)");`
92		92
93	`if(lendat > 0 ) {`	93	`if(lendat > 0 ) {`
94	`if (lendat > 1 && (nr * nc) % lendat != 0) {`	94	`if (lendat > 1 && (nr * nc) % lendat != 0) {`
95	`if (((lendat > nr) && (lendat / nr) * nr != lendat) \|\|`	95	`if (((lendat > nr) && (lendat / nr) * nr != lendat) \|\|`
96	`((lendat < nr) && (nr / lendat) * lendat != nr))`	96	`((lendat < nr) && (nr / lendat) * lendat != nr))`
97	`warning("data length [%d] is not a sub-multiple or multiple of the number of rows [%d] in matrix", lendat, nr);`	97	`warning("data length [%d] is not a sub-multiple or multiple of the number of rows [%d] in matrix", lendat, nr);`
98	`else if (((lendat > nc) && (lendat / nc) * nc != lendat) \|\|`	98	`else if (((lendat > nc) && (lendat / nc) * nc != lendat) \|\|`
99	`((lendat < nc) && (nc / lendat) * lendat != nc))`	99	`((lendat < nc) && (nc / lendat) * lendat != nc))`
100	`warning("data length [%d] is not a sub-multiple or multiple of the number of columns [%d] in matrix", lendat, nc);`	100	`warning("data length [%d] is not a sub-multiple or multiple of the number of columns [%d] in matrix", lendat, nc);`
101	`}`	101	`}`
102	`else if ((lendat > 1) && (nr * nc == 0)){`	102	`else if ((lendat > 1) && (nr * nc == 0)){`
103	`warning("data length exceeds size of matrix");`	103	`warning("data length exceeds size of matrix");`
104	`}`	104	`}`
105	`}`	105	`}`
106		106
107	`if ((double)nr * (double)nc > INT_MAX)`	107	`if ((double)nr * (double)nc > INT_MAX)`
108	`error("matrix: too many elements specified");`	108	`error("matrix: too many elements specified");`
109		109
110	`PROTECT(snr = allocMatrix(TYPEOF(vals), nr, nc));`	110	`PROTECT(snr = allocMatrix(TYPEOF(vals), nr, nc));`
111	`if(lendat) {`	111	`if(lendat) {`
112	`if (isVector(vals))`	112	`if (isVector(vals))`
113	`copyMatrix(snr, vals, byrow);`	113	`copyMatrix(snr, vals, byrow);`
114	`else`	114	`else`
115	`copyListMatrix(snr, vals, byrow);`	115	`copyListMatrix(snr, vals, byrow);`
116	`} else if (isVector(vals)) { /* fill with NAs */`	116	`} else if (isVector(vals)) { /* fill with NAs */`
117	`int i, j;`	117	`int i, j;`
118	`switch(TYPEOF(vals)) {`	118	`switch(TYPEOF(vals)) {`
119	`case STRSXP:`	119	`case STRSXP:`
120	`for (i = 0; i < nr; i++)`	120	`for (i = 0; i < nr; i++)`
121	`for (j = 0; j < nc; j++)`	121	`for (j = 0; j < nc; j++)`
122	`SET_STRING_ELT(snr, i + j * nr, NA_STRING);`	122	`SET_STRING_ELT(snr, i + j * nr, NA_STRING);`
123	`break;`	123	`break;`
124	`case LGLSXP:`	124	`case LGLSXP:`
125	`for (i = 0; i < nr; i++)`	125	`for (i = 0; i < nr; i++)`
126	`for (j = 0; j < nc; j++)`	126	`for (j = 0; j < nc; j++)`
127	`LOGICAL(snr)[i + j * nr] = NA_LOGICAL;`	127	`LOGICAL(snr)[i + j * nr] = NA_LOGICAL;`
128	`break;`	128	`break;`
129	`case INTSXP:`	129	`case INTSXP:`
130	`for (i = 0; i < nr; i++)`	130	`for (i = 0; i < nr; i++)`
131	`for (j = 0; j < nc; j++)`	131	`for (j = 0; j < nc; j++)`
132	`INTEGER(snr)[i + j * nr] = NA_INTEGER;`	132	`INTEGER(snr)[i + j * nr] = NA_INTEGER;`
133	`break;`	133	`break;`
134	`case REALSXP:`	134	`case REALSXP:`
135	`for (i = 0; i < nr; i++)`	135	`for (i = 0; i < nr; i++)`
136	`for (j = 0; j < nc; j++)`	136	`for (j = 0; j < nc; j++)`
137	`REAL(snr)[i + j * nr] = NA_REAL;`	137	`REAL(snr)[i + j * nr] = NA_REAL;`
138	`break;`	138	`break;`
139	`case CPLXSXP:`	139	`case CPLXSXP:`
140	`{`	140	`{`
141	`Rcomplex na_cmplx;`	141	`Rcomplex na_cmplx;`
142	`na_cmplx.r = NA_REAL;`	142	`na_cmplx.r = NA_REAL;`
143	`na_cmplx.i = 0;`	143	`na_cmplx.i = 0;`
144	`for (i = 0; i < nr; i++)`	144	`for (i = 0; i < nr; i++)`
145	`for (j = 0; j < nc; j++)`	145	`for (j = 0; j < nc; j++)`
146	`COMPLEX(snr)[i + j * nr] = na_cmplx;`	146	`COMPLEX(snr)[i + j * nr] = na_cmplx;`
147	`}`	147	`}`
148	`break;`	148	`break;`
-		149	`case RAWSXP:`
-		150	`for (i = 0; i < nr; i++)`
-		151	`for (j = 0; j < nc; j++)`
-		152	`RAW(snr)[i + j * nr] = 0;`
-		153	`break;`
149	`}`	154	`}`
150	`}`	155	`}`
151	`UNPROTECT(1);`	156	`UNPROTECT(1);`
152	`return snr;`	157	`return snr;`
153	`}`	158	`}`
154		159
155		160
156	`SEXP allocMatrix(SEXPTYPE mode, int nrow, int ncol)`	161	`SEXP allocMatrix(SEXPTYPE mode, int nrow, int ncol)`
157	`{`	162	`{`
158	`SEXP s, t;`	163	`SEXP s, t;`
159	`int n;`	164	`int n;`
160		165
161	`if (nrow < 0 \|\| ncol < 0)`	166	`if (nrow < 0 \|\| ncol < 0)`
162	`error("negative extents to matrix");`	167	`error("negative extents to matrix");`
163	`if ((double)nrow * (double)ncol > INT_MAX)`	168	`if ((double)nrow * (double)ncol > INT_MAX)`
164	`error("allocMatrix: too many elements specified");`	169	`error("allocMatrix: too many elements specified");`
165	`n = nrow * ncol;`	170	`n = nrow * ncol;`
166	`PROTECT(s = allocVector(mode, n));`	171	`PROTECT(s = allocVector(mode, n));`
167	`PROTECT(t = allocVector(INTSXP, 2));`	172	`PROTECT(t = allocVector(INTSXP, 2));`
168	`INTEGER(t)[0] = nrow;`	173	`INTEGER(t)[0] = nrow;`
169	`INTEGER(t)[1] = ncol;`	174	`INTEGER(t)[1] = ncol;`
170	`setAttrib(s, R_DimSymbol, t);`	175	`setAttrib(s, R_DimSymbol, t);`
171	`UNPROTECT(2);`	176	`UNPROTECT(2);`
172	`return s;`	177	`return s;`
173	`}`	178	`}`
174		179
175		180
176	`SEXP allocArray(SEXPTYPE mode, SEXP dims)`	181	`SEXP allocArray(SEXPTYPE mode, SEXP dims)`
177	`{`	182	`{`
178	`SEXP array;`	183	`SEXP array;`
179	`int i, n;`	184	`int i, n;`
180	`double dn;`	185	`double dn;`
181		186
182	`dn = n = 1;`	187	`dn = n = 1;`
183	`for (i = 0; i < LENGTH(dims); i++) {`	188	`for (i = 0; i < LENGTH(dims); i++) {`
184	`dn *= INTEGER(dims)[i];`	189	`dn *= INTEGER(dims)[i];`
185	`if(dn > INT_MAX)`	190	`if(dn > INT_MAX)`
186	`error("allocArray: too many elements specified by dims");`	191	`error("allocArray: too many elements specified by dims");`
187	`n *= INTEGER(dims)[i];`	192	`n *= INTEGER(dims)[i];`
188	`}`	193	`}`
189		194
190	`PROTECT(dims = duplicate(dims));`	195	`PROTECT(dims = duplicate(dims));`
191	`PROTECT(array = allocVector(mode, n));`	196	`PROTECT(array = allocVector(mode, n));`
192	`setAttrib(array, R_DimSymbol, dims);`	197	`setAttrib(array, R_DimSymbol, dims);`
193	`UNPROTECT(2);`	198	`UNPROTECT(2);`
194	`return array;`	199	`return array;`
195	`}`	200	`}`
196		201
197	`/* DropDims strips away redundant dimensioning information. */`	202	`/* DropDims strips away redundant dimensioning information. */`
198	`/* If there is an appropriate dimnames attribute the correct */`	203	`/* If there is an appropriate dimnames attribute the correct */`
199	`/* element is extracted and attached to the vector as a names */`	204	`/* element is extracted and attached to the vector as a names */`
200	`/* attribute. Note that this function mutates x. */`	205	`/* attribute. Note that this function mutates x. */`
201	`/* Duplication should occur before this is called. */`	206	`/* Duplication should occur before this is called. */`
202		207
203	`SEXP DropDims(SEXP x)`	208	`SEXP DropDims(SEXP x)`
204	`{`	209	`{`
205	`SEXP q, dims, dimnames, newnames = R_NilValue;`	210	`SEXP q, dims, dimnames, newnames = R_NilValue;`
206	`int i, n, ndims;`	211	`int i, n, ndims;`
207		212
208	`PROTECT(x);`	213	`PROTECT(x);`
209	`dims = getAttrib(x, R_DimSymbol);`	214	`dims = getAttrib(x, R_DimSymbol);`
210	`dimnames = getAttrib(x, R_DimNamesSymbol);`	215	`dimnames = getAttrib(x, R_DimNamesSymbol);`
211		216
212	`/* Check that dropping will actually do something. */`	217	`/* Check that dropping will actually do something. */`
213	`/* (1) Check that there is a "dim" attribute. */`	218	`/* (1) Check that there is a "dim" attribute. */`
214		219
215	`if (dims == R_NilValue) {`	220	`if (dims == R_NilValue) {`
216	`UNPROTECT(1);`	221	`UNPROTECT(1);`
217	`return x;`	222	`return x;`
218	`}`	223	`}`
219	`ndims = LENGTH(dims);`	224	`ndims = LENGTH(dims);`
220		225
221	`/* (2) Check whether there are redundant extents */`	226	`/* (2) Check whether there are redundant extents */`
222	`n = 0;`	227	`n = 0;`
223	`for (i = 0; i < ndims; i++)`	228	`for (i = 0; i < ndims; i++)`
224	`if (INTEGER(dims)[i] != 1) n++;`	229	`if (INTEGER(dims)[i] != 1) n++;`
225	`if (n == ndims) {`	230	`if (n == ndims) {`
226	`UNPROTECT(1);`	231	`UNPROTECT(1);`
227	`return x;`	232	`return x;`
228	`}`	233	`}`
229		234
230	`if (n <= 1) {`	235	`if (n <= 1) {`
231	`/* We have reduced to a vector result. */`	236	`/* We have reduced to a vector result. */`
232	`if (dimnames != R_NilValue) {`	237	`if (dimnames != R_NilValue) {`
233	`n = length(dims);`	238	`n = length(dims);`
234	`if (TYPEOF(dimnames) == VECSXP) {`	239	`if (TYPEOF(dimnames) == VECSXP) {`
235	`for (i = 0; i < n; i++) {`	240	`for (i = 0; i < n; i++) {`
236	`if (INTEGER(dims)[i] != 1) {`	241	`if (INTEGER(dims)[i] != 1) {`
237	`newnames = VECTOR_ELT(dimnames, i);`	242	`newnames = VECTOR_ELT(dimnames, i);`
238	`break;`	243	`break;`
239	`}`	244	`}`
240	`}`	245	`}`
241	`}`	246	`}`
242	`else {`	247	`else {`
243	`q = dimnames;`	248	`q = dimnames;`
244	`for (i = 0; i < n; i++) {`	249	`for (i = 0; i < n; i++) {`
245	`if (INTEGER(dims)[i] != 1) {`	250	`if (INTEGER(dims)[i] != 1) {`
246	`newnames = CAR(q);`	251	`newnames = CAR(q);`
247	`break;`	252	`break;`
248	`}`	253	`}`
249	`q = CDR(q);`	254	`q = CDR(q);`
250	`}`	255	`}`
251	`}`	256	`}`
252	`}`	257	`}`
253	`PROTECT(newnames);`	258	`PROTECT(newnames);`
254	`setAttrib(x, R_DimNamesSymbol, R_NilValue);`	259	`setAttrib(x, R_DimNamesSymbol, R_NilValue);`
255	`setAttrib(x, R_DimSymbol, R_NilValue);`	260	`setAttrib(x, R_DimSymbol, R_NilValue);`
256	`setAttrib(x, R_NamesSymbol, newnames);`	261	`setAttrib(x, R_NamesSymbol, newnames);`
257	`UNPROTECT(1);`	262	`UNPROTECT(1);`
258	`}`	263	`}`
259	`else {`	264	`else {`
260	`/* We have a lower dimensional array. */`	265	`/* We have a lower dimensional array. */`
261	`SEXP newdims, dnn, newnamesnames = R_NilValue;`	266	`SEXP newdims, dnn, newnamesnames = R_NilValue;`
262	`dnn = getAttrib(dimnames, R_NamesSymbol);`	267	`dnn = getAttrib(dimnames, R_NamesSymbol);`
263	`PROTECT(newdims = allocVector(INTSXP, n));`	268	`PROTECT(newdims = allocVector(INTSXP, n));`
264	`for (i = 0, n = 0; i < ndims; i++)`	269	`for (i = 0, n = 0; i < ndims; i++)`
265	`if (INTEGER(dims)[i] != 1)`	270	`if (INTEGER(dims)[i] != 1)`
266	`INTEGER(newdims)[n++] = INTEGER(dims)[i];`	271	`INTEGER(newdims)[n++] = INTEGER(dims)[i];`
267	`if (!isNull(dimnames)) {`	272	`if (!isNull(dimnames)) {`
268	`int havenames = 0;`	273	`int havenames = 0;`
269	`for (i = 0; i < ndims; i++)`	274	`for (i = 0; i < ndims; i++)`
270	`if (INTEGER(dims)[i] != 1 &&`	275	`if (INTEGER(dims)[i] != 1 &&`
271	`VECTOR_ELT(dimnames, i) != R_NilValue)`	276	`VECTOR_ELT(dimnames, i) != R_NilValue)`
272	`havenames = 1;`	277	`havenames = 1;`
273	`if (havenames) {`	278	`if (havenames) {`
274	`PROTECT(newnames = allocVector(VECSXP, n));`	279	`PROTECT(newnames = allocVector(VECSXP, n));`
275	`PROTECT(newnamesnames = allocVector(STRSXP, n));`	280	`PROTECT(newnamesnames = allocVector(STRSXP, n));`
276	`for (i = 0, n = 0; i < ndims; i++) {`	281	`for (i = 0, n = 0; i < ndims; i++) {`
277	`if (INTEGER(dims)[i] != 1) {`	282	`if (INTEGER(dims)[i] != 1) {`
278	`if(!isNull(dnn))`	283	`if(!isNull(dnn))`
279	`SET_STRING_ELT(newnamesnames, n,`	284	`SET_STRING_ELT(newnamesnames, n,`
280	`STRING_ELT(dnn, i));`	285	`STRING_ELT(dnn, i));`
281	`SET_VECTOR_ELT(newnames, n++, VECTOR_ELT(dimnames, i));`	286	`SET_VECTOR_ELT(newnames, n++, VECTOR_ELT(dimnames, i));`
282	`}`	287	`}`
283	`}`	288	`}`
284	`}`	289	`}`
285	`else dimnames = R_NilValue;`	290	`else dimnames = R_NilValue;`
286	`}`	291	`}`
287	`PROTECT(dimnames);`	292	`PROTECT(dimnames);`
288	`setAttrib(x, R_DimNamesSymbol, R_NilValue);`	293	`setAttrib(x, R_DimNamesSymbol, R_NilValue);`
289	`setAttrib(x, R_DimSymbol, newdims);`	294	`setAttrib(x, R_DimSymbol, newdims);`
290	`if (dimnames != R_NilValue)`	295	`if (dimnames != R_NilValue)`
291	`{`	296	`{`
292	`if(!isNull(dnn))`	297	`if(!isNull(dnn))`
293	`setAttrib(newnames, R_NamesSymbol, newnamesnames);`	298	`setAttrib(newnames, R_NamesSymbol, newnamesnames);`
294	`setAttrib(x, R_DimNamesSymbol, newnames);`	299	`setAttrib(x, R_DimNamesSymbol, newnames);`
295	`UNPROTECT(2);`	300	`UNPROTECT(2);`
296	`}`	301	`}`
297	`UNPROTECT(2);`	302	`UNPROTECT(2);`
298	`}`	303	`}`
299	`UNPROTECT(1);`	304	`UNPROTECT(1);`
300	`return x;`	305	`return x;`
301	`}`	306	`}`
302		307
303	`SEXP do_drop(SEXP call, SEXP op, SEXP args, SEXP rho)`	308	`SEXP do_drop(SEXP call, SEXP op, SEXP args, SEXP rho)`
304	`{`	309	`{`
305	`SEXP x, xdims;`	310	`SEXP x, xdims;`
306	`int i, n, shorten;`	311	`int i, n, shorten;`
307		312
308	`checkArity(op, args);`	313	`checkArity(op, args);`
309	`x = CAR(args);`	314	`x = CAR(args);`
310	`if ((xdims = getAttrib(x, R_DimSymbol)) != R_NilValue) {`	315	`if ((xdims = getAttrib(x, R_DimSymbol)) != R_NilValue) {`
311	`n = LENGTH(xdims);`	316	`n = LENGTH(xdims);`
312	`shorten = 0;`	317	`shorten = 0;`
313	`for (i = 0; i < n; i++)`	318	`for (i = 0; i < n; i++)`
314	`if (INTEGER(xdims)[i] == 1) shorten = 1;`	319	`if (INTEGER(xdims)[i] == 1) shorten = 1;`
315	`if (shorten) {`	320	`if (shorten) {`
316	`if (NAMED(x)) x = duplicate(x);`	321	`if (NAMED(x)) x = duplicate(x);`
317	`x = DropDims(x);`	322	`x = DropDims(x);`
318	`}`	323	`}`
319	`}`	324	`}`
320	`return x;`	325	`return x;`
321	`}`	326	`}`
322		327
323	`/* Length of Primitive Objects */`	328	`/* Length of Primitive Objects */`
324		329
325	`SEXP do_length(SEXP call, SEXP op, SEXP args, SEXP rho)`	330	`SEXP do_length(SEXP call, SEXP op, SEXP args, SEXP rho)`
326	`{`	331	`{`
327	`SEXP ans;`	332	`SEXP ans;`
328	`R_len_t len;`	333	`R_len_t len;`
329		334
330	`if (length(args) != 1)`	335	`if (length(args) != 1)`
331	`error("incorrect number of args to length");`	336	`error("incorrect number of args to length");`
332		337
333	`if( isObject(CAR(args)) && DispatchOrEval(call, op, "length", args,`	338	`if( isObject(CAR(args)) && DispatchOrEval(call, op, "length", args,`
334	`rho, &ans, 0, 1))`	339	`rho, &ans, 0, 1))`
335	`return(ans);`	340	`return(ans);`
336		341
337	`ans = allocVector(INTSXP, 1);`	342	`ans = allocVector(INTSXP, 1);`
338	`len = length(CAR(args));`	343	`len = length(CAR(args));`
339	`INTEGER(ans)[0] = (len < INT_MAX) ? len : NA_INTEGER;`	344	`INTEGER(ans)[0] = (len < INT_MAX) ? len : NA_INTEGER;`
340	`return ans;`	345	`return ans;`
341	`}`	346	`}`
342		347
343		348
344	`SEXP do_rowscols(SEXP call, SEXP op, SEXP args, SEXP rho)`	349	`SEXP do_rowscols(SEXP call, SEXP op, SEXP args, SEXP rho)`
345	`{`	350	`{`
346	`SEXP ans;`	351	`SEXP ans;`
347	`int i, j, nr, nc;`	352	`int i, j, nr, nc;`
348		353
349	`if (length(args) != 1)`	354	`if (length(args) != 1)`
350	`error("incorrect number of args to row/col");`	355	`error("incorrect number of args to row/col");`
351	`if (!isMatrix(CAR(args)))`	356	`if (!isMatrix(CAR(args)))`
352	`error("a matrix is required as arg to row/col");`	357	`error("a matrix is required as arg to row/col");`
353		358
354	`nr = nrows(CAR(args));`	359	`nr = nrows(CAR(args));`
355	`nc = ncols(CAR(args));`	360	`nc = ncols(CAR(args));`
356		361
357	`ans = allocMatrix(INTSXP, nr, nc);`	362	`ans = allocMatrix(INTSXP, nr, nc);`
358		363
359	`switch (PRIMVAL(op)) {`	364	`switch (PRIMVAL(op)) {`
360	`case 1:`	365	`case 1:`
361	`for (i = 0; i < nr; i++)`	366	`for (i = 0; i < nr; i++)`
362	`for (j = 0; j < nc; j++)`	367	`for (j = 0; j < nc; j++)`
363	`INTEGER(ans)[i + j * nr] = i + 1;`	368	`INTEGER(ans)[i + j * nr] = i + 1;`
364	`break;`	369	`break;`
365	`case 2:`	370	`case 2:`
366	`for (i = 0; i < nr; i++)`	371	`for (i = 0; i < nr; i++)`
367	`for (j = 0; j < nc; j++)`	372	`for (j = 0; j < nc; j++)`
368	`INTEGER(ans)[i + j * nr] = j + 1;`	373	`INTEGER(ans)[i + j * nr] = j + 1;`
369	`break;`	374	`break;`
370	`}`	375	`}`
371	`return ans;`	376	`return ans;`
372	`}`	377	`}`
373		378
374	`static void matprod(double *x, int nrx, int ncx,`	379	`static void matprod(double *x, int nrx, int ncx,`
375	`double y, int nry, int ncy, double z)`	380	`double y, int nry, int ncy, double z)`
376	`{`	381	`{`
377	`char transa = "N", transb = "N";`	382	`char transa = "N", transb = "N";`
378	`int i, j, k;`	383	`int i, j, k;`
379	`double one = 1.0, zero = 0.0, sum;`	384	`double one = 1.0, zero = 0.0, sum;`
380	`Rboolean have_na = FALSE;`	385	`Rboolean have_na = FALSE;`
381		386
382	`if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {`	387	`if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {`
383	`/* Don't trust the BLAS to handle NA/NaNs correctly: PR#4582`	388	`/* Don't trust the BLAS to handle NA/NaNs correctly: PR#4582`
384	`* The test is only O(n) here`	389	`* The test is only O(n) here`
385	`*/`	390	`*/`
386	`for (i = 0; i < nrx*ncx; i++)`	391	`for (i = 0; i < nrx*ncx; i++)`
387	`if (ISNAN(x[i])) {have_na = TRUE; break;}`	392	`if (ISNAN(x[i])) {have_na = TRUE; break;}`
388	`if (!have_na)`	393	`if (!have_na)`
389	`for (i = 0; i < nry*ncy; i++)`	394	`for (i = 0; i < nry*ncy; i++)`
390	`if (ISNAN(y[i])) {have_na = TRUE; break;}`	395	`if (ISNAN(y[i])) {have_na = TRUE; break;}`
391	`if (have_na) {`	396	`if (have_na) {`
392	`for (i = 0; i < nrx; i++)`	397	`for (i = 0; i < nrx; i++)`
393	`for (k = 0; k < ncy; k++) {`	398	`for (k = 0; k < ncy; k++) {`
394	`sum = 0.0;`	399	`sum = 0.0;`
395	`for (j = 0; j < ncx; j++)`	400	`for (j = 0; j < ncx; j++)`
396	`sum += x[i + j * nrx] * y[j + k * nry];`	401	`sum += x[i + j * nrx] * y[j + k * nry];`
397	`z[i + k * nrx] = sum;`	402	`z[i + k * nrx] = sum;`
398	`}`	403	`}`
399	`} else`	404	`} else`
400	`F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,`	405	`F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,`
401	`x, &nrx, y, &nry, &zero, z, &nrx);`	406	`x, &nrx, y, &nry, &zero, z, &nrx);`
402	`} else /* zero-extent operations should return zeroes */`	407	`} else /* zero-extent operations should return zeroes */`
403	`for(i = 0; i < nrx*ncy; i++) z[i] = 0;`	408	`for(i = 0; i < nrx*ncy; i++) z[i] = 0;`
404	`}`	409	`}`
405		410
406	`#ifdef HAVE_DOUBLE_COMPLEX`	411	`#ifdef HAVE_DOUBLE_COMPLEX`
407	`/* ZGEMM - perform one of the matrix-matrix operations */`	412	`/* ZGEMM - perform one of the matrix-matrix operations */`
408	`/* C := alphaop( A )op( B ) + betaC /`	413	`/* C := alphaop( A )op( B ) + betaC /`
409	`extern void`	414	`extern void`
410	`F77_NAME(zgemm)(const char transa, const char transb, const int *m,`	415	`F77_NAME(zgemm)(const char transa, const char transb, const int *m,`
411	`const int n, const int k, const Rcomplex *alpha,`	416	`const int n, const int k, const Rcomplex *alpha,`
412	`const Rcomplex a, const int lda,`	417	`const Rcomplex a, const int lda,`
413	`const Rcomplex b, const int ldb,`	418	`const Rcomplex b, const int ldb,`
414	`const Rcomplex beta, Rcomplex c, const int *ldc);`	419	`const Rcomplex beta, Rcomplex c, const int *ldc);`
415	`#endif`	420	`#endif`
416		421
417	`static void cmatprod(Rcomplex *x, int nrx, int ncx,`	422	`static void cmatprod(Rcomplex *x, int nrx, int ncx,`
418	`Rcomplex y, int nry, int ncy, Rcomplex z)`	423	`Rcomplex y, int nry, int ncy, Rcomplex z)`
419	`{`	424	`{`
420	`#ifdef HAVE_DOUBLE_COMPLEX`	425	`#ifdef HAVE_DOUBLE_COMPLEX`
421	`char transa = "N", transb = "N";`	426	`char transa = "N", transb = "N";`
422	`int i;`	427	`int i;`
423	`Rcomplex one, zero;`	428	`Rcomplex one, zero;`
424		429
425	`one.r = 1.0; one.i = zero.r = zero.i = 0.0;`	430	`one.r = 1.0; one.i = zero.r = zero.i = 0.0;`
426	`if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {`	431	`if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {`
427	`F77_CALL(zgemm)(transa, transb, &nrx, &ncy, &ncx, &one,`	432	`F77_CALL(zgemm)(transa, transb, &nrx, &ncy, &ncx, &one,`
428	`x, &nrx, y, &nry, &zero, z, &nrx);`	433	`x, &nrx, y, &nry, &zero, z, &nrx);`
429	`} else { /* zero-extent operations should return zeroes */`	434	`} else { /* zero-extent operations should return zeroes */`
430	`for(i = 0; i < nrx*ncy; i++) z[i].r = z[i].i = 0;`	435	`for(i = 0; i < nrx*ncy; i++) z[i].r = z[i].i = 0;`
431	`}`	436	`}`
432	`#else`	437	`#else`
433	`int i, j, k;`	438	`int i, j, k;`
434	`double xij_r, xij_i, yjk_r, yjk_i, sum_i, sum_r;`	439	`double xij_r, xij_i, yjk_r, yjk_i, sum_i, sum_r;`
435		440
436	`for (i = 0; i < nrx; i++)`	441	`for (i = 0; i < nrx; i++)`
437	`for (k = 0; k < ncy; k++) {`	442	`for (k = 0; k < ncy; k++) {`
438	`z[i + k * nrx].r = NA_REAL;`	443	`z[i + k * nrx].r = NA_REAL;`
439	`z[i + k * nrx].i = NA_REAL;`	444	`z[i + k * nrx].i = NA_REAL;`
440	`sum_r = 0.0;`	445	`sum_r = 0.0;`
441	`sum_i = 0.0;`	446	`sum_i = 0.0;`
442	`for (j = 0; j < ncx; j++) {`	447	`for (j = 0; j < ncx; j++) {`
443	`xij_r = x[i + j * nrx].r;`	448	`xij_r = x[i + j * nrx].r;`
444	`xij_i = x[i + j * nrx].i;`	449	`xij_i = x[i + j * nrx].i;`
445	`yjk_r = y[j + k * nry].r;`	450	`yjk_r = y[j + k * nry].r;`
446	`yjk_i = y[j + k * nry].i;`	451	`yjk_i = y[j + k * nry].i;`
447	`if (ISNAN(xij_r) \|\| ISNAN(xij_i)`	452	`if (ISNAN(xij_r) \|\| ISNAN(xij_i)`
448	`\|\| ISNAN(yjk_r) \|\| ISNAN(yjk_i))`	453	`\|\| ISNAN(yjk_r) \|\| ISNAN(yjk_i))`
449	`goto next_ik;`	454	`goto next_ik;`
450	`sum_r += (xij_r * yjk_r - xij_i * yjk_i);`	455	`sum_r += (xij_r * yjk_r - xij_i * yjk_i);`
451	`sum_i += (xij_r * yjk_i + xij_i * yjk_r);`	456	`sum_i += (xij_r * yjk_i + xij_i * yjk_r);`
452	`}`	457	`}`
453	`z[i + k * nrx].r = sum_r;`	458	`z[i + k * nrx].r = sum_r;`
454	`z[i + k * nrx].i = sum_i;`	459	`z[i + k * nrx].i = sum_i;`
455	`next_ik:`	460	`next_ik:`
456	`;`	461	`;`
457	`}`	462	`}`
458	`#endif`	463	`#endif`
459	`}`	464	`}`
460		465
461	`static void symcrossprod(double x, int nr, int nc, double z)`	466	`static void symcrossprod(double x, int nr, int nc, double z)`
462	`{`	467	`{`
463	`char trans = "T", uplo = "U";`	468	`char trans = "T", uplo = "U";`
464	`double one = 1.0, zero = 0.0;`	469	`double one = 1.0, zero = 0.0;`
465	`int i, j;`	470	`int i, j;`
466	`if (nr > 0 && nc > 0) {`	471	`if (nr > 0 && nc > 0) {`
467	`F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc);`	472	`F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc);`
468	`for (i = 1; i < nc; i++)`	473	`for (i = 1; i < nc; i++)`
469	`for (j = 0; j < i; j++) z[i + nc j] = z[j + nc i];`	474	`for (j = 0; j < i; j++) z[i + nc j] = z[j + nc i];`
470	`} else { /* zero-extent operations should return zeroes */`	475	`} else { /* zero-extent operations should return zeroes */`
471	`for(i = 0; i < nc*nc; i++) z[i] = 0;`	476	`for(i = 0; i < nc*nc; i++) z[i] = 0;`
472	`}`	477	`}`
473		478
474	`}`	479	`}`
475		480
476	`static void crossprod(double *x, int nrx, int ncx,`	481	`static void crossprod(double *x, int nrx, int ncx,`
477	`double y, int nry, int ncy, double z)`	482	`double y, int nry, int ncy, double z)`
478	`{`	483	`{`
479	`char transa = "T", transb = "N";`	484	`char transa = "T", transb = "N";`
480	`double one = 1.0, zero = 0.0;`	485	`double one = 1.0, zero = 0.0;`
481	`if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {`	486	`if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {`
482	`F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one,`	487	`F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one,`
483	`x, &nrx, y, &nry, &zero, z, &ncx);`	488	`x, &nrx, y, &nry, &zero, z, &ncx);`
484	`}`	489	`}`
485	`else { /* zero-extent operations should return zeroes */`	490	`else { /* zero-extent operations should return zeroes */`
486	`int i;`	491	`int i;`
487	`for(i = 0; i < ncx*ncy; i++) z[i] = 0;`	492	`for(i = 0; i < ncx*ncy; i++) z[i] = 0;`
488	`}`	493	`}`
489	`}`	494	`}`
490		495
491	`static void ccrossprod(Rcomplex *x, int nrx, int ncx,`	496	`static void ccrossprod(Rcomplex *x, int nrx, int ncx,`
492	`Rcomplex y, int nry, int ncy, Rcomplex z)`	497	`Rcomplex y, int nry, int ncy, Rcomplex z)`
493	`{`	498	`{`
494	`char transa = "T", transb = "N";`	499	`char transa = "T", transb = "N";`
495	`Rcomplex one, zero;`	500	`Rcomplex one, zero;`

The R Project SVN R

(root)/branches/djm-source/src/main/array.c @ 34412 – Rev 29864 → 30691