The R Project SVN R

Rev

Rev 78127 | Rev 80325 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 78127 Rev 80272
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  R : A Computer Language for Statistical Data Analysis
 2 
 *  R : A Computer Language for Statistical Data Analysis
3 
 *  Copyright (C) 1998-2020   The R Core Team
 3 
 *  Copyright (C) 1998-2021   The R Core Team
4 
 *  Copyright (C) 2002-2015   The R Foundation
 4 
 *  Copyright (C) 2002-2015   The R Foundation
5 
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
 5 
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
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
Line 61... Line 61...
61 
	return VECTOR_ELT(dimnames, 1);
 61 
	return VECTOR_ELT(dimnames, 1);
62 
    else
 62 
    else
63 
	return R_NilValue;
 63 
	return R_NilValue;
64 
}
 64 
}
65 
 
 65 
 
- 
 
 66 
// .Internal(matrix(data, nrow, ncol, byrow, dimnames, missing(nrow), missing(ncol)))
66 
SEXP attribute_hidden do_matrix(SEXP call, SEXP op, SEXP args, SEXP rho)
 67 
SEXP attribute_hidden do_matrix(SEXP call, SEXP op, SEXP args, SEXP rho)
67 
{
 68 
{
68 
    SEXP vals, ans, snr, snc, dimnames;
 69 
    SEXP vals, ans, snr, snc, dimnames;
69 
    int nr = 1, nc = 1, byrow, miss_nr, miss_nc;
 70 
    int nr = 1, nc = 1, byrow, miss_nr, miss_nc;
70 
    R_xlen_t lendat;
 71 
    R_xlen_t lendat;
Line 114... Line 115...
114 
    }
 115 
    }
115 
    if (miss_nr && miss_nc) {
 116 
    if (miss_nr && miss_nc) {
116 
	if (lendat > INT_MAX) error("data is too long");
 117 
	if (lendat > INT_MAX) error("data is too long");
117 
	nr = (int) lendat;
 118 
	nr = (int) lendat;
118 
    } else if (miss_nr) {
 119 
    } else if (miss_nr) {
119 
	if (lendat > (double) nc * INT_MAX) error("data is too long");
 120 
	if (lendat > (double) nc * INT_MAX) error("data is too long"); // incl lendat > nc == 0
120 
	// avoid division by zero
 - 
 
121 
	if (nc == 0) {
 - 
 
122 
	    if (lendat) error(_("nc = 0 for non-null data"));
 121 
	if (nc == 0) // as lendat <= nc, have lendat == 0
123 
	    else nr = 0;
 122 
	    nr = 0;
124 
	} else
 123 
	else
125 
	    nr = (int) ceil((double) lendat / (double) nc);
 124 
	    nr = (int) ceil((double) lendat / (double) nc);
126 
    } else if (miss_nc) {
 125 
    } else if (miss_nc) {
127 
	if (lendat > (double) nr * INT_MAX) error("data is too long");
 126 
	if (lendat > (double) nr * INT_MAX) error("data is too long"); // incl lendat > nr == 0
128 
	// avoid division by zero
 - 
 
129 
	if (nr == 0) {
 127 
	if (nr == 0) // then lendat == 0
130 
	    if (lendat) error(_("nr = 0 for non-null data"));
 - 
 
131 
	    else nc = 0;
 128 
	    nc = 0;
132 
	} else
 129 
	else
133 
	    nc = (int) ceil((double) lendat / (double) nr);
 130 
	    nc = (int) ceil((double) lendat / (double) nr);
134 
    }
 131 
    }
135 
 
 132 
 
136 
    if(lendat > 0) {
 133 
    if (lendat > 0) {
137 
	R_xlen_t nrc = (R_xlen_t) nr * nc;
 134 
	R_xlen_t nrc = (R_xlen_t) nr * nc;
138 
	if (lendat > 1 && nrc % lendat != 0) {
 135 
	if (lendat > 1 && (nrc % lendat) != 0) { // ==> nrc > 0
139 
	    if (((lendat > nr) && (lendat / nr) * nr != lendat) ||
 136 
	    if (((lendat > nr) && (lendat / nr) * nr != lendat) ||
140 
		((lendat < nr) && (nr / lendat) * lendat != nr))
 137 
		((lendat < nr) && (nr / lendat) * lendat != nr))
141 
		warning(_("data length [%d] is not a sub-multiple or multiple of the number of rows [%d]"), lendat, nr);
 138 
		warning(_("data length [%d] is not a sub-multiple or multiple of the number of rows [%d]"), lendat, nr);
142 
	    else if (((lendat > nc) && (lendat / nc) * nc != lendat) ||
 139 
	    else if (((lendat > nc) && (lendat / nc) * nc != lendat) ||
143 
		     ((lendat < nc) && (nc / lendat) * lendat != nc))
 140 
		     ((lendat < nc) && (nc / lendat) * lendat != nc))
144 
		warning(_("data length [%d] is not a sub-multiple or multiple of the number of columns [%d]"), lendat, nc);
 141 
		warning(_("data length [%d] is not a sub-multiple or multiple of the number of columns [%d]"), lendat, nc);
- 
 
 142 
	    else if (nrc != lendat)
- 
 
 143 
		warning(_("data length differs from size of matrix: [%d != %d x %d]"), lendat, nr, nc);
145 
	}
 144 
	}
146 
	else if ((lendat > 1) && (nrc == 0)){
 145 
	else if (lendat > 1 && nrc == 0) // for now *not* warning for e.g., matrix(NA, 0, 4)
147 
	    warning(_("data length exceeds size of matrix"));
 146 
	    warning(_("non-empty data for zero-extent matrix"));
148 
	}
 - 
 
149 
    }
 147 
    }
150 
 
 148 
 
151 
#ifndef LONG_VECTOR_SUPPORT
 149 
#ifndef LONG_VECTOR_SUPPORT
152 
    if ((double)nr * (double)nc > INT_MAX)
 150 
    if ((double)nr * (double)nc > INT_MAX)
153 
	error(_("too many elements specified"));
 151 
	error(_("too many elements specified"));
Line 287... Line 285...
287 
/* attribute.  Note that this function mutates x. */
 285 
/* attribute.  Note that this function mutates x. */
288 
/* Duplication should occur before this is called. */
 286 
/* Duplication should occur before this is called. */
289 
 
 287 
 
290 
SEXP DropDims(SEXP x)
 288 
SEXP DropDims(SEXP x)
291 
{
 289 
{
292 
    SEXP dims, dimnames, newnames = R_NilValue;
 - 
 
293 
    int i, n, ndims;
 - 
 
294 
 
 - 
 
295 
    PROTECT(x);
 290 
    PROTECT(x);
296 
    dims = getAttrib(x, R_DimSymbol);
 291 
    SEXP dims = getAttrib(x, R_DimSymbol);
297 
 
 292 
 
298 
    /* Check that dropping will actually do something. */
 293 
    /* Check that dropping will actually do something. */
299 
    /* (1) Check that there is a "dim" attribute. */
 294 
    /* (1) Check that there is a "dim" attribute. */
300 
 
 295 
 
301 
    if (dims == R_NilValue) {
 296 
    if (dims == R_NilValue) {
302 
	UNPROTECT(1); /* x */
 297 
	UNPROTECT(1); /* x */
303 
	return x;
 298 
	return x;
304 
    }
 299 
    }
- 
 
 300 
 
305 
    ndims = LENGTH(dims);
 301 
    int ndims = LENGTH(dims);
306 
    int *dim = INTEGER(dims); // used several times
 302 
    int *dim = INTEGER(dims); // used several times
307 
 
 303 
 
308 
    /* (2) Check whether there are redundant extents */
 304 
    /* (2) Check whether there are redundant extents */
309 
    n = 0;
 305 
    int i, n = 0;
310 
    for (i = 0; i < ndims; i++)
 306 
    for (i = 0; i < ndims; i++)
311 
	if (dim[i] != 1) n++;
 307 
	if (dim[i] != 1) n++;
312 
    if (n == ndims) {
 308 
    if (n == ndims) {
313 
	UNPROTECT(1); /* x */
 309 
	UNPROTECT(1); /* x */
314 
	return x;
 310 
	return x;
315 
    }
 311 
    }
316 
 
 312 
 
317 
    PROTECT(dimnames = getAttrib(x, R_DimNamesSymbol));
 313 
    SEXP dimnames = PROTECT(getAttrib(x, R_DimNamesSymbol)),
- 
 
 314 
	newnames = R_NilValue;
318 
    if (n <= 1) {
 315 
    if (n <= 1) {
319 
	/* We have reduced to a vector result.
 316 
	/* We have reduced to a vector result.
320 
	   If that has length one, it is ambiguous which dimnames to use,
 317 
	   If that has length one, it is ambiguous which dimnames to use,
321 
	   so use it if there is only one (as from R 2.7.0).
 318 
	   so use it if there is only one (as from R 2.7.0).
322 
	 */
 319 
	 */
Line 809... Line 806...
809 
	/* Instead of xY, compute (xY)^T == (Y^T)(x^T)
 806 
	/* Instead of xY, compute (xY)^T == (Y^T)(x^T)
810 
	   The result is a vector, so transposing its content is no-op */
 807 
	   The result is a vector, so transposing its content is no-op */
811 
	F77_CALL(dgemv)(transT, &nry, &ncy, &one, y,
 808 
	F77_CALL(dgemv)(transT, &nry, &ncy, &one, y,
812 
			&nry, x, &ione, &zero, z, &ione FCONE);
 809 
			&nry, x, &ione, &zero, z, &ione FCONE);
813 
    else /* matrix-matrix or outer product */
 810 
    else /* matrix-matrix or outer product */
814 
	F77_CALL(dgemm)(transN, transN, &nrx, &ncy, &ncx, &one, x,
811 
	F77_CALL(dgemm)(transN, transN, &nrx, &ncy, &ncx, &one, x,
815 
			&nrx, y, &nry, &zero, z, &nrx FCONE FCONE);
 812 
			&nrx, y, &nry, &zero, z, &nrx FCONE FCONE);
816 
}
 813 
}
817 
 
 814 
 
818 
static void internal_cmatprod(Rcomplex *x, int nrx, int ncx,
 815 
static void internal_cmatprod(Rcomplex *x, int nrx, int ncx,
819 
                              Rcomplex *y, int nry, int ncy, Rcomplex *z)
 816 
                              Rcomplex *y, int nry, int ncy, Rcomplex *z)