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

 /*
  *  R : A Computer Language for Statistical Data Analysis
  *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
- *  Copyright (C) 1998-2003   The R Development Core Team.
+ *  Copyright (C) 1998-2001   The R Development Core Team
+ *  Copyright (C) 2002--2003  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.
+ *
- *  You should have received a copy of the GNU General Public License
+ *  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
- *  along with this program; if not, write to the Free Software
+ *  writing to the Free Software Foundation, Inc., 59 Temple Place,
- *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ *  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 */
-/* FIXME : This is nonsense.  When the "dimnames" attribute is */
-/* grabbed off an array it is always adjusted to be a vector. */
+/* "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 if (TYPEOF(dimnames) == LISTSXP)
-	return CAR(dimnames);
     else
 	return R_NilValue;
+}
 SEXP GetColNames(SEXP dimnames)
+{
     if (TYPEOF(dimnames) == VECSXP)
 	return VECTOR_ELT(dimnames, 1);
-    else if (TYPEOF(dimnames) == LISTSXP)
-	return CADR(dimnames);
     else
 	return R_NilValue;
+}
 SEXP do_matrix(SEXP call, SEXP op, SEXP args, SEXP rho)
+}
 #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);
     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];
+    }
+}
 static void crossprod(double *x, int nrx, int ncx,
+	    ;
 #endif
+	}
 #endif
+}
+/* "%*%" (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;
     xdims = getAttrib(x, R_DimSymbol);
     ydims = getAttrib(y, R_DimSymbol);
     ldx = length(xdims);
     ldy = length(ydims);
-    if (ldx != 2 && ldy != 2) {
+    if (ldx != 2 && ldy != 2) {		/* x and y non-matrices */
 	if (PRIMVAL(op) == 0) {
 	    nrx = 1;
 	    ncx = LENGTH(x);
+	}
 	else {
 	    ncx = 1;
+	}
 	nry = LENGTH(y);
 	ncy = 1;
+    }
-    else if (ldx != 2) {
+    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 */
+	    if (LENGTH(x) == nry) {	/* x as row vector */
 		nrx = 1;
 		ncx = LENGTH(x);
+	    }
-	    else if (nry == 1) { /* x as col vector */
+	    else if (nry == 1) {	/* x as col vector */
 		nrx = LENGTH(x);
 		ncx = 1;
+	    }
+	}
 	else {
-	    if (LENGTH(x) == nry) {
+	    if (LENGTH(x) == nry) {	/* x is a row vector */
 		nrx = LENGTH(x);
 		ncx = 1;
+	    }
+	}
+    }
-    else if (ldy != 2) {
+    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 */
+	    if (LENGTH(y) == ncx) {	/* y as col vector */
 		nry = LENGTH(y);
 		ncy = 1;
+	    }
-	    else if (ncx == 1){ /* y as row vector */
+	    else if (ncx == 1) {	/* y as row vector */
 		nry = 1;
 		ncy = LENGTH(y);
+	    }
+	}
 	else {
-	    if (LENGTH(y) == nrx) {
+	    if (LENGTH(y) == nrx) {	/* y is a row vector */
 		nry = LENGTH(y);
 		ncy = 1;
+	    }
+	}
+    }
-    else {
+    else {				/* x and y matrices */
 	nrx = INTEGER(xdims)[0];
 	ncx = INTEGER(xdims)[1];
 	nry = INTEGER(ydims)[0];
 	ncy = INTEGER(ydims)[1];
+    }
     else
 	mode = REALSXP;
     SETCAR(args, coerceVector(CAR(args), mode));
     SETCADR(args, coerceVector(CADR(args), mode));
-    if (PRIMVAL(op) == 0) {
+    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) {
 	    setAttrib(dimnames, R_NamesSymbol, dimnamesnames);
 	    setAttrib(ans, R_DimNamesSymbol, dimnames);
 	    UNPROTECT(2);
+	}
+    }
-    else {
+    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
 #endif
 		crossprod(REAL(CAR(args)), nrx, ncx,
 			  REAL(CADR(args)), nry, ncy, REAL(ans));
+	}
 	PROTECT(xdims = getAttrib(CAR(args), R_DimNamesSymbol));
-	if(sym)
+	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);
+		    dnx = getAttrib(xdims, R_NamesSymbol);
-		SET_VECTOR_ELT(dimnames, 0, VECTOR_ELT(xdims, 1));
+		    SET_VECTOR_ELT(dimnames, 0, VECTOR_ELT(xdims, 1));
-		if(!isNull(dnx))
+		    if(!isNull(dnx))
-		    SET_STRING_ELT(dimnamesnames, 0, STRING_ELT(dnx, 1));
+			SET_STRING_ELT(dimnamesnames, 0, STRING_ELT(dnx, 1));
+		}
+	    }
 	    if (ydims != R_NilValue) {
+                if (ldy == 2) {
-		dny = getAttrib(ydims, R_NamesSymbol);
+		    dny = getAttrib(ydims, R_NamesSymbol);
-		SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 1));
+		    SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 1));
-		if(!isNull(dny))
+		    if(!isNull(dny))
-		    SET_STRING_ELT(dimnamesnames, 1, STRING_ELT(dny, 1));
+			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);
+	}
 	    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);

The R Project SVN R

(root)/branches/R-grid/src/main/array.c – Rev 26266 → 26606