| /* |
| * R : A Computer Language for Statistical Data Analysis |
| * Copyright (C) 1998-2018 The R Core Team |
| * Copyright (C) 2002-2015 The R Foundation |
| * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka |
| * |
| * 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 |
| * along with this program; if not, a copy is available at |
| * https://www.R-project.org/Licenses/ |
| */ |
| |
| #ifdef HAVE_CONFIG_H |
| #include <config.h> |
| #endif |
| |
| #include <Defn.h> |
| #include <Internal.h> |
| #include <Rmath.h> |
| #include <R_ext/RS.h> /* for Calloc/Free */ |
| #include <R_ext/Applic.h> /* for dgemm */ |
| #include <R_ext/Itermacros.h> |
| |
| #include "duplicate.h" |
| |
| #include <complex.h> |
| #include "Rcomplex.h" /* toC99 */ |
| |
| /* "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.'' |
| |
| They are used in bind.c and subset.c, and advertised in Rinternals.h |
| */ |
| 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 attribute_hidden do_matrix(SEXP call, SEXP op, SEXP args, SEXP rho) |
| { |
| SEXP vals, ans, snr, snc, dimnames; |
| int nr = 1, nc = 1, byrow, miss_nr, miss_nc; |
| R_xlen_t lendat; |
| |
| checkArity(op, args); |
| vals = CAR(args); args = CDR(args); |
| switch(TYPEOF(vals)) { |
| case LGLSXP: |
| case INTSXP: |
| case REALSXP: |
| case CPLXSXP: |
| case STRSXP: |
| case RAWSXP: |
| case EXPRSXP: |
| case VECSXP: |
| break; |
| default: |
| error(_("'data' must be of a vector type, was '%s'"), |
| type2char(TYPEOF(vals))); |
| } |
| lendat = XLENGTH(vals); |
| snr = CAR(args); args = CDR(args); |
| snc = CAR(args); args = CDR(args); |
| byrow = asLogical(CAR(args)); args = CDR(args); |
| if (byrow == NA_INTEGER) |
| error(_("invalid '%s' argument"), "byrow"); |
| dimnames = CAR(args); |
| args = CDR(args); |
| miss_nr = asLogical(CAR(args)); args = CDR(args); |
| miss_nc = asLogical(CAR(args)); |
| |
| if (!miss_nr) { |
| if (!isNumeric(snr)) error(_("non-numeric matrix extent")); |
| nr = asInteger(snr); |
| if (nr == NA_INTEGER) |
| error(_("invalid 'nrow' value (too large or NA)")); |
| if (nr < 0) |
| error(_("invalid 'nrow' value (< 0)")); |
| } |
| if (!miss_nc) { |
| if (!isNumeric(snc)) error(_("non-numeric matrix extent")); |
| nc = asInteger(snc); |
| if (nc == NA_INTEGER) |
| error(_("invalid 'ncol' value (too large or NA)")); |
| if (nc < 0) |
| error(_("invalid 'ncol' value (< 0)")); |
| } |
| if (miss_nr && miss_nc) { |
| if (lendat > INT_MAX) error("data is too long"); |
| nr = (int) lendat; |
| } else if (miss_nr) { |
| if (lendat > (double) nc * INT_MAX) error("data is too long"); |
| // avoid division by zero |
| if (nc == 0) { |
| if (lendat) error(_("nc = 0 for non-null data")); |
| else nr = 0; |
| } else |
| nr = (int) ceil((double) lendat / (double) nc); |
| } else if (miss_nc) { |
| if (lendat > (double) nr * INT_MAX) error("data is too long"); |
| // avoid division by zero |
| if (nr == 0) { |
| if (lendat) error(_("nr = 0 for non-null data")); |
| else nc = 0; |
| } else |
| nc = (int) ceil((double) lendat / (double) nr); |
| } |
| |
| if(lendat > 0) { |
| R_xlen_t nrc = (R_xlen_t) nr * nc; |
| if (lendat > 1 && nrc % 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]"), 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]"), lendat, nc); |
| } |
| else if ((lendat > 1) && (nrc == 0)){ |
| warning(_("data length exceeds size of matrix")); |
| } |
| } |
| |
| #ifndef LONG_VECTOR_SUPPORT |
| if ((double)nr * (double)nc > INT_MAX) |
| error(_("too many elements specified")); |
| #endif |
| |
| PROTECT(ans = allocMatrix(TYPEOF(vals), nr, nc)); |
| if(lendat) { |
| if (isVector(vals)) |
| copyMatrix(ans, vals, byrow); |
| else |
| copyListMatrix(ans, vals, byrow); |
| } else if (isVector(vals)) { /* fill with NAs */ |
| R_xlen_t N = (R_xlen_t) nr * nc, i; |
| switch(TYPEOF(vals)) { |
| case STRSXP: |
| for (i = 0; i < N; i++) |
| SET_STRING_ELT(ans, i, NA_STRING); |
| break; |
| case LGLSXP: |
| for (i = 0; i < N; i++) |
| LOGICAL(ans)[i] = NA_LOGICAL; |
| break; |
| case INTSXP: |
| for (i = 0; i < N; i++) |
| INTEGER(ans)[i] = NA_INTEGER; |
| break; |
| case REALSXP: |
| for (i = 0; i < N; i++) |
| REAL(ans)[i] = NA_REAL; |
| break; |
| case CPLXSXP: |
| { |
| Rcomplex na_cmplx; |
| na_cmplx.r = NA_REAL; |
| na_cmplx.i = 0; |
| for (i = 0; i < N; i++) |
| COMPLEX(ans)[i] = na_cmplx; |
| } |
| break; |
| case RAWSXP: |
| if (N) memset(RAW(ans), 0, N); |
| break; |
| default: |
| /* don't fill with anything */ |
| ; |
| } |
| } |
| if(!isNull(dimnames) && length(dimnames) > 0) |
| ans = dimnamesgets(ans, dimnames); |
| UNPROTECT(1); |
| return ans; |
| } |
| |
| |
| SEXP allocMatrix(SEXPTYPE mode, int nrow, int ncol) |
| { |
| SEXP s, t; |
| R_xlen_t n; |
| |
| if (nrow < 0 || ncol < 0) |
| error(_("negative extents to matrix")); |
| #ifndef LONG_VECTOR_SUPPORT |
| if ((double)nrow * (double)ncol > INT_MAX) |
| error(_("allocMatrix: too many elements specified")); |
| #endif |
| n = ((R_xlen_t) 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; |
| } |
| |
| /** |
| * Allocate a 3-dimensional array |
| * |
| * @param mode The R mode (e.g. INTSXP) |
| * @param nrow number of rows |
| * @param ncol number of columns |
| * @param nface number of faces |
| * |
| * @return A 3-dimensional array of the indicated dimensions and mode |
| */ |
| SEXP alloc3DArray(SEXPTYPE mode, int nrow, int ncol, int nface) |
| { |
| SEXP s, t; |
| R_xlen_t n; |
| |
| if (nrow < 0 || ncol < 0 || nface < 0) |
| error(_("negative extents to 3D array")); |
| #ifndef LONG_VECTOR_SUPPORT |
| if ((double)nrow * (double)ncol * (double)nface > INT_MAX) |
| error(_("'alloc3DArray': too many elements specified")); |
| #endif |
| n = ((R_xlen_t) nrow) * ncol * nface; |
| PROTECT(s = allocVector(mode, n)); |
| PROTECT(t = allocVector(INTSXP, 3)); |
| INTEGER(t)[0] = nrow; |
| INTEGER(t)[1] = ncol; |
| INTEGER(t)[2] = nface; |
| setAttrib(s, R_DimSymbol, t); |
| UNPROTECT(2); |
| return s; |
| } |
| |
| |
| SEXP allocArray(SEXPTYPE mode, SEXP dims) |
| { |
| SEXP array; |
| int i; |
| R_xlen_t n = 1; |
| #ifndef LONG_VECTOR_SUPPORT |
| double dn = 1; |
| #endif |
| |
| for (i = 0; i < LENGTH(dims); i++) { |
| #ifndef LONG_VECTOR_SUPPORT |
| dn *= INTEGER(dims)[i]; |
| if(dn > INT_MAX) |
| error(_("'allocArray': too many elements specified by 'dims'")); |
| #endif |
| 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 dims, dimnames, newnames = R_NilValue; |
| int i, n, ndims; |
| |
| PROTECT(x); |
| dims = getAttrib(x, R_DimSymbol); |
| |
| /* Check that dropping will actually do something. */ |
| /* (1) Check that there is a "dim" attribute. */ |
| |
| if (dims == R_NilValue) { |
| UNPROTECT(1); /* x */ |
| return x; |
| } |
| ndims = LENGTH(dims); |
| int *dim = INTEGER(dims); // used several times |
| |
| /* (2) Check whether there are redundant extents */ |
| n = 0; |
| for (i = 0; i < ndims; i++) |
| if (dim[i] != 1) n++; |
| if (n == ndims) { |
| UNPROTECT(1); /* x */ |
| return x; |
| } |
| |
| PROTECT(dimnames = getAttrib(x, R_DimNamesSymbol)); |
| if (n <= 1) { |
| /* We have reduced to a vector result. |
| If that has length one, it is ambiguous which dimnames to use, |
| so use it if there is only one (as from R 2.7.0). |
| */ |
| if (dimnames != R_NilValue) { |
| if(XLENGTH(x) != 1) { |
| for (i = 0; i < LENGTH(dims); i++) { |
| if (dim[i] != 1) { |
| newnames = VECTOR_ELT(dimnames, i); |
| break; |
| } |
| } |
| } else { /* drop all dims: keep names if unambiguous */ |
| int cnt; |
| for(i = 0, cnt = 0; i < LENGTH(dims); i++) |
| if(VECTOR_ELT(dimnames, i) != R_NilValue) cnt++; |
| if(cnt == 1) |
| for (i = 0; i < LENGTH(dims); i++) { |
| newnames = VECTOR_ELT(dimnames, i); |
| if(newnames != R_NilValue) break; |
| } |
| } |
| } |
| PROTECT(newnames); |
| setAttrib(x, R_DimNamesSymbol, R_NilValue); |
| setAttrib(x, R_DimSymbol, R_NilValue); |
| setAttrib(x, R_NamesSymbol, newnames); |
| /* FIXME: the following is desirable, but pointless as long as |
| subset.c & others have a contrary version that leaves the |
| S4 class in, incorrectly, in the case of vectors. JMC |
| 3/3/09 */ |
| /* if(IS_S4_OBJECT(x)) {/\* no longer valid subclass of array or |
| matrix *\/ */ |
| /* setAttrib(x, R_ClassSymbol, R_NilValue); */ |
| /* UNSET_S4_OBJECT(x); */ |
| /* } */ |
| UNPROTECT(1); /* newnames */ |
| } else { |
| // We have a lower dimensional array, and n == length(newdims) |
| SEXP newdims, dnn, newnamesnames = R_NilValue; |
| PROTECT(dnn = getAttrib(dimnames, R_NamesSymbol)); |
| PROTECT(newdims = allocVector(INTSXP, n)); |
| for (i = 0, n = 0; i < ndims; i++) |
| if (dim[i] != 1) |
| INTEGER(newdims)[n++] = dim[i]; |
| if(!isNull(getAttrib(dims, R_NamesSymbol))) { |
| SEXP nms_d = getAttrib(dims, R_NamesSymbol), |
| new_nms = PROTECT(allocVector(STRSXP, n)); |
| for (i = 0, n = 0; i < ndims; i++) |
| if (dim[i] != 1) |
| SET_STRING_ELT(new_nms, n++, STRING_ELT(nms_d, i)); |
| setAttrib(newdims, R_NamesSymbol, new_nms); |
| UNPROTECT(1); |
| } |
| Rboolean havenames = FALSE; |
| if (!isNull(dimnames)) { |
| for (i = 0; i < ndims; i++) |
| if (dim[i] != 1 && |
| VECTOR_ELT(dimnames, i) != R_NilValue) |
| havenames = TRUE; |
| if (havenames) { |
| PROTECT(newnames = allocVector(VECSXP, n)); |
| PROTECT(newnamesnames = allocVector(STRSXP, n)); |
| for (i = 0, n = 0; i < ndims; i++) { |
| if (dim[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; |
| } |
| setAttrib(x, R_DimNamesSymbol, R_NilValue); |
| setAttrib(x, R_DimSymbol, newdims); |
| if (havenames) |
| { |
| if(!isNull(dnn)) |
| setAttrib(newnames, R_NamesSymbol, newnamesnames); |
| setAttrib(x, R_DimNamesSymbol, newnames); |
| UNPROTECT(2); /* newnamesnames, newnames */ |
| } |
| UNPROTECT(2); /* newdims, dnn */ |
| } |
| UNPROTECT(2); /* dimnames, x */ |
| return x; |
| } |
| |
| SEXP attribute_hidden 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 (MAYBE_REFERENCED(x)) x = R_duplicate_attr(x); |
| x = DropDims(x); |
| } |
| } |
| return x; |
| } |
| |
| /* Length of Primitive Objects */ |
| |
| SEXP attribute_hidden do_length(SEXP call, SEXP op, SEXP args, SEXP rho) |
| { |
| checkArity(op, args); |
| check1arg(args, call, "x"); |
| |
| SEXP x = CAR(args), ans; |
| |
| if (isObject(x) && |
| DispatchOrEval(call, op, "length", args, rho, &ans, 0, 1)) { |
| if (length(ans) == 1 && TYPEOF(ans) == REALSXP) { |
| double d = REAL(ans)[0]; |
| if (R_FINITE(d) && d >= 0. && d <= INT_MAX && floor(d) == d) { |
| PROTECT(ans); |
| ans = coerceVector(ans, INTSXP); |
| UNPROTECT(1); |
| return(ans); |
| } |
| } |
| return(ans); |
| } |
| |
| |
| #ifdef LONG_VECTOR_SUPPORT |
| // or use IS_LONG_VEC |
| R_xlen_t len = xlength(x); |
| if (len > INT_MAX) return ScalarReal((double) len); |
| #endif |
| return ScalarInteger(length(x)); |
| } |
| |
| R_len_t attribute_hidden dispatch_length(SEXP x, SEXP call, SEXP rho) { |
| R_xlen_t len = dispatch_xlength(x, call, rho); |
| #ifdef LONG_VECTOR_SUPPORT |
| if (len > INT_MAX) return R_BadLongVector(x, __FILE__, __LINE__); |
| #endif |
| return (R_len_t) len; |
| } |
| |
| R_xlen_t attribute_hidden dispatch_xlength(SEXP x, SEXP call, SEXP rho) { |
| static SEXP length_op = NULL; |
| if (isObject(x)) { |
| SEXP len, args; |
| if (length_op == NULL) |
| length_op = R_Primitive("length"); |
| PROTECT(args = list1(x)); |
| if (DispatchOrEval(call, length_op, "length", args, rho, &len, 0, 1)) { |
| UNPROTECT(1); |
| return (R_xlen_t) |
| (TYPEOF(len) == REALSXP ? REAL(len)[0] : asInteger(len)); |
| } |
| UNPROTECT(1); |
| } |
| return(xlength(x)); |
| } |
| |
| // auxiliary for do_lengths_*(), i.e., R's lengths() |
| static R_xlen_t getElementLength(SEXP x, R_xlen_t i, SEXP call, SEXP rho) { |
| SEXP x_elt; |
| R_xlen_t ans; |
| |
| PROTECT(x_elt = dispatch_subset2(x, i, call, rho)); |
| ans = dispatch_xlength(x_elt, call, rho); |
| UNPROTECT(1); /* x_elt */ |
| return ans; |
| } |
| |
| #ifdef LONG_VECTOR_SUPPORT |
| static SEXP do_lengths_long(SEXP x, SEXP call, SEXP rho) |
| { |
| SEXP ans; |
| R_xlen_t x_len, i; |
| double *ans_elt; |
| |
| x_len = dispatch_xlength(x, call, rho); |
| PROTECT(ans = allocVector(REALSXP, x_len)); |
| for (i = 0, ans_elt = REAL(ans); i < x_len; i++, ans_elt++) |
| *ans_elt = (double) getElementLength(x, i, call, rho); |
| UNPROTECT(1); |
| return ans; |
| } |
| #endif |
| |
| SEXP attribute_hidden do_lengths(SEXP call, SEXP op, SEXP args, SEXP rho) |
| { |
| checkArity(op, args); |
| SEXP x = CAR(args), ans; |
| R_xlen_t x_len, i; |
| int *ans_elt; |
| int useNames = asLogical(CADR(args)); |
| if (useNames == NA_LOGICAL) |
| error(_("invalid '%s' value"), "use.names"); |
| |
| if (DispatchOrEval(call, op, "lengths", args, rho, &ans, 0, 1)) |
| return(ans); |
| |
| Rboolean isList = isVectorList(x) || isS4(x); |
| if(!isList) switch(TYPEOF(x)) { |
| case NILSXP: |
| case CHARSXP: |
| case LGLSXP: |
| case INTSXP: |
| case REALSXP: |
| case CPLXSXP: |
| case STRSXP: |
| case RAWSXP: |
| break; |
| default: |
| error(_("'%s' must be a list or atomic vector"), "x"); |
| } |
| x_len = dispatch_xlength(x, call, rho); |
| PROTECT(ans = allocVector(INTSXP, x_len)); |
| if(isList) { |
| for (i = 0, ans_elt = INTEGER(ans); i < x_len; i++, ans_elt++) { |
| R_xlen_t x_elt_len = getElementLength(x, i, call, rho); |
| #ifdef LONG_VECTOR_SUPPORT |
| if (x_elt_len > INT_MAX) { |
| ans = do_lengths_long(x, call, rho); |
| UNPROTECT(1); |
| PROTECT(ans); |
| break; |
| } |
| #endif |
| *ans_elt = (int)x_elt_len; |
| } |
| } else { // atomic: every element has length 1 |
| for (i = 0, ans_elt = INTEGER(ans); i < x_len; i++, ans_elt++) |
| *ans_elt = 1; |
| } |
| SEXP dim = getAttrib(x, R_DimSymbol); |
| if(!isNull(dim)) { |
| setAttrib(ans, R_DimSymbol, dim); |
| } |
| if(useNames) { |
| SEXP names = getAttrib(x, R_NamesSymbol); |
| if(!isNull(names)) setAttrib(ans, R_NamesSymbol, names); |
| SEXP dimnames = getAttrib(x, R_DimNamesSymbol); |
| if(!isNull(dimnames)) setAttrib(ans, R_DimNamesSymbol, dimnames); |
| } |
| UNPROTECT(1); |
| return ans; |
| } |
| |
| SEXP attribute_hidden do_rowscols(SEXP call, SEXP op, SEXP args, SEXP rho) |
| { |
| checkArity(op, args); |
| SEXP dim = CAR(args); |
| int nprot = 0; |
| if (!isInteger(dim)) { |
| PROTECT(dim = coerceVector(dim, INTSXP)); nprot++; |
| } |
| if (LENGTH(dim) != 2) |
| error(_("a matrix-like object is required as argument to '%s'"), |
| (PRIMVAL(op) == 2) ? "col" : "row"); |
| |
| int nr = INTEGER(dim)[0], |
| nc = INTEGER(dim)[1]; |
| if(nprot) UNPROTECT(nprot); |
| |
| SEXP ans = allocMatrix(INTSXP, nr, nc); |
| |
| R_xlen_t NR = nr; |
| switch (PRIMVAL(op)) { |
| case 1: // row() & .row() |
| for (int i = 0; i < nr; i++) |
| for (int j = 0; j < nc; j++) |
| INTEGER(ans)[i + j * NR] = i + 1; |
| break; |
| case 2: // col() & .col() |
| for (int i = 0; i < nr; i++) |
| for (int j = 0; j < nc; j++) |
| INTEGER(ans)[i + j * NR] = j + 1; |
| break; |
| } |
| return ans; |
| } |
| |
| /* |
| Whenever vector x contains NaN or Inf (or -Inf), the function returns TRUE. |
| It can be imprecise: it can return TRUE in other cases as well. |
| |
| A precise version of the function could be implemented as |
| |
| for (R_xlen_t i = 0; i < n; i++) |
| if (!R_FINITE(x[i])) return TRUE; |
| return FALSE; |
| |
| The present version is imprecise, but faster. |
| */ |
| static Rboolean mayHaveNaNOrInf(double *x, R_xlen_t n) |
| { |
| if ((n&1) != 0 && !R_FINITE(x[0])) |
| return TRUE; |
| for (R_xlen_t i = n&1; i < n; i += 2) |
| /* A precise version could use this condition: |
| * |
| * !R_FINITE(x[i]+x[i+1]) && (!R_FINITE(x[i]) || !R_FINITE(x[i+1])) |
| * |
| * The present imprecise version has been found to be faster |
| * with GCC and ICC in the common case when the sum of the two |
| * values is always finite. |
| * |
| * The present version is imprecise because the sum of two very |
| * large finite values (e.g. 1e308) may be infinite. |
| */ |
| if (!R_FINITE(x[i]+x[i+1])) |
| return TRUE; |
| return FALSE; |
| } |
| |
| /* |
| This is an experimental version that has been observed to run fast on some |
| SIMD hardware with GCC and ICC. |
| Note that the OpenMP reduction assumes associativity of addition, which is |
| safe here, because the result is only used for an imprecise test for |
| the presence of NaN and Inf values. |
| */ |
| static Rboolean mayHaveNaNOrInf_simd(double *x, R_xlen_t n) |
| { |
| double s = 0; |
| /* SIMD reduction is supported since OpenMP 4.0. The value of _OPENMP is |
| unreliable in some compilers, so we depend on HAVE_OPENMP_SIMDRED, |
| which is normally set by configure based on a test. */ |
| /* _OPENMP >= 201307 */ |
| #if defined(_OPENMP) && HAVE_OPENMP_SIMDRED |
| #pragma omp simd reduction(+:s) |
| #endif |
| for (R_xlen_t i = 0; i < n; i++) |
| s += x[i]; |
| return !R_FINITE(s); |
| } |
| |
| static Rboolean cmayHaveNaNOrInf(Rcomplex *x, R_xlen_t n) |
| { |
| /* With HAVE_FORTRAN_DOUBLE_COMPLEX set, it should be clear that |
| Rcomplex has no padding, so we could probably use mayHaveNaNOrInf, |
| but better safe than sorry... */ |
| if ((n&1) != 0 && (!R_FINITE(x[0].r) || !R_FINITE(x[0].i))) |
| return TRUE; |
| for (R_xlen_t i = n&1; i < n; i += 2) |
| if (!R_FINITE(x[i].r+x[i].i+x[i+1].r+x[i+1].i)) |
| return TRUE; |
| return FALSE; |
| } |
| |
| /* experimental version for SIMD hardware (see also mayHaveNaNOrInf_simd) */ |
| static Rboolean cmayHaveNaNOrInf_simd(Rcomplex *x, R_xlen_t n) |
| { |
| double s = 0; |
| /* _OPENMP >= 201307 - see mayHaveNaNOrInf_simd */ |
| #if defined(_OPENMP) && HAVE_OPENMP_SIMDRED |
| #pragma omp simd reduction(+:s) |
| #endif |
| for (R_xlen_t i = 0; i < n; i++) { |
| s += x[i].r; |
| s += x[i].i; |
| } |
| return !R_FINITE(s); |
| } |
| |
| static void internal_matprod(double *x, int nrx, int ncx, |
| double *y, int nry, int ncy, double *z) |
| { |
| LDOUBLE sum; |
| #define MATPROD_BODY \ |
| R_xlen_t NRX = nrx, NRY = nry; \ |
| for (int i = 0; i < nrx; i++) \ |
| for (int k = 0; k < ncy; k++) { \ |
| sum = 0.0; \ |
| for (int j = 0; j < ncx; j++) \ |
| sum += x[i + j * NRX] * y[j + k * NRY]; \ |
| z[i + k * NRX] = (double) sum; \ |
| } |
| MATPROD_BODY; |
| } |
| |
| static void simple_matprod(double *x, int nrx, int ncx, |
| double *y, int nry, int ncy, double *z) |
| { |
| double sum; |
| MATPROD_BODY; |
| } |
| |
| static void internal_crossprod(double *x, int nrx, int ncx, |
| double *y, int nry, int ncy, double *z) |
| { |
| LDOUBLE sum; |
| #define CROSSPROD_BODY \ |
| R_xlen_t NRX = nrx, NRY = nry, NCX = ncx; \ |
| for (int i = 0; i < ncx; i++) \ |
| for (int k = 0; k < ncy; k++) { \ |
| sum = 0.0; \ |
| for (int j = 0; j < nrx; j++) \ |
| sum += x[j + i * NRX] * y[j + k * NRY]; \ |
| z[i + k * NCX] = (double) sum; \ |
| } |
| CROSSPROD_BODY; |
| } |
| |
| static void simple_crossprod(double *x, int nrx, int ncx, |
| double *y, int nry, int ncy, double *z) |
| { |
| double sum; |
| CROSSPROD_BODY; |
| } |
| |
| static void internal_tcrossprod(double *x, int nrx, int ncx, |
| double *y, int nry, int ncy, double *z) |
| { |
| LDOUBLE sum; |
| #define TCROSSPROD_BODY \ |
| R_xlen_t NRX = nrx, NRY = nry; \ |
| for (int i = 0; i < nrx; i++) \ |
| for (int k = 0; k < nry; k++) { \ |
| sum = 0.0; \ |
| for (int j = 0; j < ncx; j++) \ |
| sum += x[i + j * NRX] * y[k + j * NRY]; \ |
| z[i + k * NRX] = (double) sum; \ |
| } |
| TCROSSPROD_BODY; |
| } |
| |
| static void simple_tcrossprod(double *x, int nrx, int ncx, |
| double *y, int nry, int ncy, double *z) |
| { |
| double sum; |
| TCROSSPROD_BODY; |
| } |
| |
| |
| static void matprod(double *x, int nrx, int ncx, |
| double *y, int nry, int ncy, double *z) |
| { |
| R_xlen_t NRX = nrx, NRY = nry; |
| if (nrx == 0 || ncx == 0 || nry == 0 || ncy == 0) { |
| /* zero-extent operations should return zeroes */ |
| for(R_xlen_t i = 0; i < NRX*ncy; i++) z[i] = 0; |
| return; |
| } |
| |
| switch(R_Matprod) { |
| case MATPROD_DEFAULT: |
| /* Don't trust the BLAS to handle NA/NaNs correctly: PR#4582 |
| * The test is only O(n) here. |
| * |
| * MKL disclaimer: "LAPACK routines assume that input matrices |
| * do not contain IEEE 754 special values such as INF or NaN values. |
| * Using these special values may cause LAPACK to return unexpected |
| * results or become unstable." |
| */ |
| if (mayHaveNaNOrInf(x, NRX*ncx) || mayHaveNaNOrInf(y, NRY*ncy)) { |
| simple_matprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| } |
| break; /* use blas */ |
| case MATPROD_INTERNAL: |
| internal_matprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| case MATPROD_BLAS: |
| break; |
| case MATPROD_DEFAULT_SIMD: |
| if (mayHaveNaNOrInf_simd(x, NRX*ncx) || |
| mayHaveNaNOrInf_simd(y, NRY*ncy)) { |
| simple_matprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| } |
| break; /* use blas */ |
| } |
| |
| char *transN = "N", *transT = "T"; |
| double one = 1.0, zero = 0.0; |
| int ione = 1; |
| |
| if (ncy == 1) /* matrix-vector or dot product */ |
| F77_CALL(dgemv)(transN, &nrx, &ncx, &one, x, |
| &nrx, y, &ione, &zero, z, &ione FCONE); |
| else if (nrx == 1) /* vector-matrix */ |
| /* Instead of xY, compute (xY)^T == (Y^T)(x^T) |
| The result is a vector, so transposing its content is no-op */ |
| F77_CALL(dgemv)(transT, &nry, &ncy, &one, y, |
| &nry, x, &ione, &zero, z, &ione FCONE); |
| else /* matrix-matrix or outer product */ |
| F77_CALL(dgemm)(transN, transN, &nrx, &ncy, &ncx, &one, x, |
| &nrx, y, &nry, &zero, z, &nrx FCONE FCONE); |
| } |
| |
| static void internal_cmatprod(Rcomplex *x, int nrx, int ncx, |
| Rcomplex *y, int nry, int ncy, Rcomplex *z) |
| { |
| LDOUBLE sum_i, sum_r; |
| #define CMATPROD_BODY \ |
| int i, j, k; \ |
| double complex xij, yjk; \ |
| R_xlen_t NRX = nrx, NRY = nry; \ |
| for (i = 0; i < nrx; i++) \ |
| for (k = 0; k < ncy; k++) { \ |
| sum_r = 0.0; \ |
| sum_i = 0.0; \ |
| for (j = 0; j < ncx; j++) { \ |
| xij = toC99(x + (i + j * NRX)); \ |
| yjk = toC99(y + (j + k * NRY)); \ |
| sum_r += creal(xij * yjk); \ |
| sum_i += cimag(xij * yjk); \ |
| } \ |
| z[i + k * NRX].r = (double) sum_r; \ |
| z[i + k * NRX].i = (double) sum_i; \ |
| } |
| CMATPROD_BODY; |
| } |
| |
| static void simple_cmatprod(Rcomplex *x, int nrx, int ncx, |
| Rcomplex *y, int nry, int ncy, Rcomplex *z) |
| { |
| double sum_i, sum_r; |
| CMATPROD_BODY; |
| } |
| |
| static void internal_ccrossprod(Rcomplex *x, int nrx, int ncx, |
| Rcomplex *y, int nry, int ncy, Rcomplex *z) |
| { |
| LDOUBLE sum_i, sum_r; |
| #define CCROSSPROD_BODY \ |
| int i, j, k; \ |
| double complex xji, yjk; \ |
| R_xlen_t NRX = nrx, NRY = nry, NCX = ncx; \ |
| for (i = 0; i < ncx; i++) \ |
| for (k = 0; k < ncy; k++) { \ |
| sum_r = 0.0; \ |
| sum_i = 0.0; \ |
| for (j = 0; j < nrx; j++) { \ |
| xji = toC99(x + (j + i * NRX)); \ |
| yjk = toC99(y + (j + k * NRY)); \ |
| sum_r += creal(xji * yjk); \ |
| sum_i += cimag(xji * yjk); \ |
| } \ |
| z[i + k * NCX].r = (double) sum_r; \ |
| z[i + k * NCX].i = (double) sum_i; \ |
| } |
| CCROSSPROD_BODY; |
| } |
| |
| static void simple_ccrossprod(Rcomplex *x, int nrx, int ncx, |
| Rcomplex *y, int nry, int ncy, Rcomplex *z) |
| { |
| double sum_i, sum_r; |
| CCROSSPROD_BODY; |
| } |
| |
| static void internal_tccrossprod(Rcomplex *x, int nrx, int ncx, |
| Rcomplex *y, int nry, int ncy, Rcomplex *z) |
| { |
| LDOUBLE sum_i, sum_r; |
| #define TCCROSSPROD_BODY \ |
| int i, j, k; \ |
| double complex xij, ykj; \ |
| R_xlen_t NRX = nrx, NRY = nry; \ |
| for (i = 0; i < nrx; i++) \ |
| for (k = 0; k < nry; k++) { \ |
| sum_r = 0.0; \ |
| sum_i = 0.0; \ |
| for (j = 0; j < ncx; j++) { \ |
| xij = toC99(x + (i + j * NRX)); \ |
| ykj = toC99(y + (k + j * NRY)); \ |
| sum_r += creal(xij * ykj); \ |
| sum_i += cimag(xij * ykj); \ |
| } \ |
| z[i + k * NRX].r = (double) sum_r; \ |
| z[i + k * NRX].i = (double) sum_i; \ |
| } |
| TCCROSSPROD_BODY; |
| } |
| |
| static void simple_tccrossprod(Rcomplex *x, int nrx, int ncx, |
| Rcomplex *y, int nry, int ncy, Rcomplex *z) |
| { |
| double sum_i, sum_r; |
| TCCROSSPROD_BODY; |
| } |
| |
| static void cmatprod(Rcomplex *x, int nrx, int ncx, |
| Rcomplex *y, int nry, int ncy, Rcomplex *z) |
| { |
| R_xlen_t NRX = nrx, NRY = nry; |
| if (nrx == 0 || ncx == 0 || nry == 0 || ncy == 0) { |
| /* zero-extent operations should return zeroes */ |
| for(R_xlen_t i = 0; i < NRX*ncy; i++) z[i].r = z[i].i = 0; |
| return; |
| } |
| |
| #ifndef HAVE_FORTRAN_DOUBLE_COMPLEX |
| if (R_Matprod == MATPROD_INTERNAL) |
| internal_cmatprod(x, nrx, ncx, y, nry, ncy, z); |
| else |
| simple_cmatprod(x, nrx, ncx, y, nry, ncy, z); |
| #else |
| switch(R_Matprod) { |
| case MATPROD_DEFAULT: |
| if (cmayHaveNaNOrInf(x, NRX*ncx) || cmayHaveNaNOrInf(y, NRY*ncy)) { |
| simple_cmatprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| } |
| break; /* use blas */ |
| case MATPROD_INTERNAL: |
| internal_cmatprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| case MATPROD_BLAS: |
| break; |
| case MATPROD_DEFAULT_SIMD: |
| if (cmayHaveNaNOrInf_simd(x, NRX*ncx) || |
| cmayHaveNaNOrInf_simd(y, NRY*ncy)) { |
| simple_cmatprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| } |
| break; /* use blas */ |
| } |
| |
| char *transa = "N", *transb = "N"; |
| Rcomplex one, zero; |
| one.r = 1.0; one.i = zero.r = zero.i = 0.0; |
| |
| F77_CALL(zgemm)(transa, transb, &nrx, &ncy, &ncx, &one, |
| x, &nrx, y, &nry, &zero, z, &nrx FCONE FCONE); |
| #endif |
| } |
| |
| static void symcrossprod(double *x, int nr, int nc, double *z) |
| { |
| R_xlen_t NR = nr, NC = nc; |
| if (nr == 0 || nc == 0) { |
| /* zero-extent operations should return zeroes */ |
| for(R_xlen_t i = 0; i < NC*NC; i++) z[i] = 0; |
| return; |
| } |
| |
| switch(R_Matprod) { |
| case MATPROD_DEFAULT: |
| /* see matprod for more details */ |
| if (mayHaveNaNOrInf(x, NR*nc)) { |
| simple_crossprod(x, nr, nc, x, nr, nc, z); |
| return; |
| } |
| break; /* use blas */ |
| case MATPROD_INTERNAL: |
| internal_crossprod(x, nr, nc, x, nr, nc, z); |
| return; |
| case MATPROD_BLAS: |
| break; |
| case MATPROD_DEFAULT_SIMD: |
| if (mayHaveNaNOrInf_simd(x, NR*nc)) { |
| simple_crossprod(x, nr, nc, x, nr, nc, z); |
| return; |
| } |
| break; /* use blas */ |
| } |
| |
| char *trans = "T", *uplo = "U"; |
| double one = 1.0, zero = 0.0; |
| |
| F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc |
| FCONE FCONE); |
| for (int i = 1; i < nc; i++) |
| for (int j = 0; j < i; j++) z[i + NC *j] = z[j + NC * i]; |
| } |
| |
| static void crossprod(double *x, int nrx, int ncx, |
| double *y, int nry, int ncy, double *z) |
| { |
| R_xlen_t NRX = nrx, NRY = nry; |
| if (nrx == 0 || ncx == 0 || nry == 0 || ncy == 0) { |
| /* zero-extent operations should return zeroes */ |
| R_xlen_t NCX = ncx; |
| for(R_xlen_t i = 0; i < NCX*ncy; i++) z[i] = 0; |
| return; |
| } |
| |
| switch(R_Matprod) { |
| case MATPROD_DEFAULT: |
| /* see matprod for more details */ |
| if (mayHaveNaNOrInf(x, NRX*ncx) || mayHaveNaNOrInf(y, NRY*ncy)) { |
| simple_crossprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| } |
| break; /* use blas */ |
| case MATPROD_INTERNAL: |
| internal_crossprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| case MATPROD_BLAS: |
| break; |
| case MATPROD_DEFAULT_SIMD: |
| if (mayHaveNaNOrInf_simd(x, NRX*ncx) || |
| mayHaveNaNOrInf_simd(y, NRY*ncy)) { |
| simple_crossprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| } |
| break; /* use blas */ |
| } |
| |
| char *transT = "T", *transN = "N"; |
| double one = 1.0, zero = 0.0; |
| int ione = 1; |
| |
| if (ncy == 1) /* matrix-vector or dot product */ |
| F77_CALL(dgemv)(transT, &nrx, &ncx, &one, x, |
| &nrx, y, &ione, &zero, z, &ione FCONE); |
| else if (ncx == 1) /* vector-matrix */ |
| /* Instead of (x^T)Y, compute ((x^T)Y)^T == (Y^T)x |
| The result is a vector, so transposing its content is no-op */ |
| F77_CALL(dgemv)(transT, &nry, &ncy, &one, y, |
| &nry, x, &ione, &zero, z, &ione FCONE); |
| else /* matrix-matrix or outer product */ |
| F77_CALL(dgemm)(transT, transN, &ncx, &ncy, &nrx, &one, |
| x, &nrx, y, &nry, &zero, z, &ncx FCONE FCONE); |
| } |
| |
| static void ccrossprod(Rcomplex *x, int nrx, int ncx, |
| Rcomplex *y, int nry, int ncy, Rcomplex *z) |
| { |
| R_xlen_t NRX = nrx, NRY = nry; |
| if (nrx == 0 || ncx == 0 || nry == 0 || ncy == 0) { |
| /* zero-extent operations should return zeroes */ |
| R_xlen_t NCX = ncx; |
| for(R_xlen_t i = 0; i < NCX*ncy; i++) z[i].r = z[i].i = 0; |
| return; |
| } |
| |
| #ifndef HAVE_FORTRAN_DOUBLE_COMPLEX |
| if (R_Matprod == MATPROD_INTERNAL) |
| internal_ccrossprod(x, nrx, ncx, y, nry, ncy, z); |
| else |
| simple_ccrossprod(x, nrx, ncx, y, nry, ncy, z); |
| #else |
| switch(R_Matprod) { |
| case MATPROD_DEFAULT: |
| if (cmayHaveNaNOrInf(x, NRX*ncx) || cmayHaveNaNOrInf(y, NRY*ncy)) { |
| simple_ccrossprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| } |
| break; /* use blas */ |
| case MATPROD_INTERNAL: |
| internal_ccrossprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| case MATPROD_BLAS: |
| break; |
| case MATPROD_DEFAULT_SIMD: |
| if (cmayHaveNaNOrInf_simd(x, NRX*ncx) || |
| cmayHaveNaNOrInf_simd(y, NRY*ncy)) { |
| simple_ccrossprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| } |
| break; /* use blas */ |
| } |
| |
| char *transa = "T", *transb = "N"; |
| Rcomplex one, zero; |
| one.r = 1.0; one.i = zero.r = zero.i = 0.0; |
| |
| F77_CALL(zgemm)(transa, transb, &ncx, &ncy, &nrx, &one, |
| x, &nrx, y, &nry, &zero, z, &ncx FCONE FCONE); |
| #endif |
| } |
| |
| static void symtcrossprod(double *x, int nr, int nc, double *z) |
| { |
| R_xlen_t NR = nr; |
| if (nr == 0 || nc == 0) { |
| /* zero-extent operations should return zeroes */ |
| for(R_xlen_t i = 0; i < NR*NR; i++) z[i] = 0; |
| return; |
| } |
| |
| switch(R_Matprod) { |
| case MATPROD_DEFAULT: |
| /* see matprod for more details */ |
| if (mayHaveNaNOrInf(x, NR*nc)) { |
| simple_tcrossprod(x, nr, nc, x, nr, nc, z); |
| return; |
| } |
| break; /* use blas */ |
| case MATPROD_INTERNAL: |
| internal_tcrossprod(x, nr, nc, x, nr, nc, z); |
| return; |
| case MATPROD_BLAS: |
| break; |
| case MATPROD_DEFAULT_SIMD: |
| if (mayHaveNaNOrInf_simd(x, NR*nc)) { |
| simple_tcrossprod(x, nr, nc, x, nr, nc, z); |
| return; |
| } |
| break; /* use blas */ |
| } |
| |
| char *trans = "N", *uplo = "U"; |
| double one = 1.0, zero = 0.0; |
| |
| F77_CALL(dsyrk)(uplo, trans, &nr, &nc, &one, x, &nr, &zero, z, &nr |
| FCONE FCONE); |
| for (int i = 1; i < nr; i++) |
| for (int j = 0; j < i; j++) z[i + nr *j] = z[j + nr * i]; |
| } |
| |
| static void tcrossprod(double *x, int nrx, int ncx, |
| double *y, int nry, int ncy, double *z) |
| { |
| R_xlen_t NRX = nrx, NRY = nry; |
| if (nrx == 0 || ncx == 0 || nry == 0 || ncy == 0) { |
| /* zero-extent operations should return zeroes */ |
| for(R_xlen_t i = 0; i < NRX*nry; i++) z[i] = 0; |
| return; |
| } |
| |
| switch(R_Matprod) { |
| case MATPROD_DEFAULT: |
| if (mayHaveNaNOrInf(x, NRX*ncx) || mayHaveNaNOrInf(y, NRY*ncy)) { |
| simple_tcrossprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| } |
| break; /* use blas */ |
| case MATPROD_INTERNAL: |
| internal_tcrossprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| case MATPROD_BLAS: |
| break; |
| case MATPROD_DEFAULT_SIMD: |
| if (mayHaveNaNOrInf_simd(x, NRX*ncx) || |
| mayHaveNaNOrInf_simd(y, NRY*ncy)) { |
| simple_tcrossprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| } |
| break; /* use blas */ |
| } |
| |
| char *transN = "N", *transT = "T"; |
| double one = 1.0, zero = 0.0; |
| int ione = 1; |
| |
| if (nry == 1) /* matrix-vector or dot product */ |
| F77_CALL(dgemv)(transN, &nrx, &ncx, &one, x, |
| &nrx, y, &ione, &zero, z, &ione FCONE); |
| else if (nrx == 1) /* vector-matrix */ |
| /* Instead of x(Y^T), compute (x(Y^T))^T == Y(x^T) |
| The result is a vector, so transposing its content is no-op */ |
| F77_CALL(dgemv)(transN, &nry, &ncy, &one, y, |
| &nry, x, &ione, &zero, z, &ione FCONE); |
| else /* matrix-matrix or outer product */ |
| F77_CALL(dgemm)(transN, transT, &nrx, &nry, &ncx, &one, |
| x, &nrx, y, &nry, &zero, z, &nrx FCONE FCONE); |
| } |
| |
| static void tccrossprod(Rcomplex *x, int nrx, int ncx, |
| Rcomplex *y, int nry, int ncy, Rcomplex *z) |
| { |
| R_xlen_t NRX = nrx, NRY = nry; |
| if (nrx == 0 || ncx == 0 || nry == 0 || ncy == 0) { |
| /* zero-extent operations should return zeroes */ |
| for(R_xlen_t i = 0; i < NRX*nry; i++) z[i].r = z[i].i = 0; |
| return; |
| } |
| |
| #ifndef HAVE_FORTRAN_DOUBLE_COMPLEX |
| if (R_Matprod == MATPROD_INTERNAL) |
| internal_tccrossprod(x, nrx, ncx, y, nry, ncy, z); |
| else |
| simple_tccrossprod(x, nrx, ncx, y, nry, ncy, z); |
| #else |
| switch(R_Matprod) { |
| case MATPROD_DEFAULT: |
| if (cmayHaveNaNOrInf(x, NRX*ncx) || cmayHaveNaNOrInf(y, NRY*ncy)) { |
| simple_tccrossprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| } |
| break; /* use blas */ |
| case MATPROD_INTERNAL: |
| internal_tccrossprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| case MATPROD_BLAS: |
| break; |
| case MATPROD_DEFAULT_SIMD: |
| if (cmayHaveNaNOrInf_simd(x, NRX*ncx) || |
| cmayHaveNaNOrInf_simd(y, NRY*ncy)) { |
| simple_tccrossprod(x, nrx, ncx, y, nry, ncy, z); |
| return; |
| } |
| break; /* use blas */ |
| } |
| |
| char *transa = "N", *transb = "T"; |
| Rcomplex one, zero; |
| one.r = 1.0; one.i = zero.r = zero.i = 0.0; |
| |
| F77_CALL(zgemm)(transa, transb, &nrx, &nry, &ncx, &one, |
| x, &nrx, y, &nry, &zero, z, &nrx FCONE FCONE); |
| #endif |
| } |
| |
| |
| /* "%*%" (op = 0), crossprod (op = 1) or tcrossprod (op = 2) */ |
| SEXP attribute_hidden 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 (PRIMVAL(op) == 0 && /* %*% is primitive, the others are .Internal() */ |
| (IS_S4_OBJECT(x) || IS_S4_OBJECT(y)) |
| && R_has_methods(op)) { |
| SEXP s, value; |
| /* Remove argument names to ensure positional matching */ |
| for(s = args; s != R_NilValue; s = CDR(s)) SET_TAG(s, R_NilValue); |
| value = R_possible_dispatch(call, op, args, rho, FALSE); |
| if (value) return value; |
| } |
| |
| checkArity(op, args); |
| sym = isNull(y); |
| if (sym && (PRIMVAL(op) > 0)) y = x; |
| if ( !(isNumeric(x) || isComplex(x)) || !(isNumeric(y) || isComplex(y)) ) |
| errorcall(call, _("requires numeric/complex 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 */ |
| // for crossprod, allow two cases: n x n ==> (1,n) x (n,1); 1 x n = (n, 1) x (1, n) |
| if (PRIMVAL(op) == 1 && LENGTH(x) == 1) { |
| nrx = ncx = nry = 1; |
| ncy = LENGTH(y); |
| } |
| else { |
| nry = LENGTH(y); |
| ncy = 1; |
| if (PRIMVAL(op) == 0) { |
| nrx = 1; |
| ncx = LENGTH(x); |
| if(ncx == 1) { // y as row vector |
| ncy = nry; |
| nry = 1; |
| } |
| } |
| else { |
| nrx = LENGTH(x); |
| ncx = 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 = nry; /* == LENGTH(x) */ |
| } |
| else if (nry == 1) { /* x as col vector */ |
| nrx = LENGTH(x); |
| ncx = 1; |
| } |
| } |
| else if (PRIMVAL(op) == 1) { /* crossprod() */ |
| if (LENGTH(x) == nry) { /* x is a col vector */ |
| nrx = nry; /* == LENGTH(x) */ |
| ncx = 1; |
| } |
| /* else if (nry == 1) ... not being too tolerant |
| to treat x as row vector, as t(x) *is* row vector */ |
| } |
| else { /* tcrossprod */ |
| if (LENGTH(x) == ncy) { /* x as row vector */ |
| nrx = 1; |
| ncx = ncy; /* == LENGTH(x) */ |
| } |
| else if (ncy == 1) { /* x as col 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 = ncx; |
| ncy = 1; |
| } |
| else if (ncx == 1) { /* y as row vector */ |
| nry = 1; |
| ncy = LENGTH(y); |
| } |
| } |
| else if (PRIMVAL(op) == 1) { /* crossprod() */ |
| if (LENGTH(y) == nrx) { /* y is a col vector */ |
| nry = nrx; |
| ncy = 1; |
| } else if (nrx == 1) { // y as row vector |
| nry = 1; |
| ncy = LENGTH(y); |
| } |
| } |
| else { // tcrossprod |
| if (nrx == 1) { // y as row vector |
| nry = 1; |
| ncy = LENGTH(y); |
| } |
| else { // y is a col 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]; |
| } |
| /* nr[ow](.) and nc[ol](.) are now defined for x and y */ |
| |
| if (PRIMVAL(op) == 0) { |
| /* primitive, so use call */ |
| if (ncx != nry) |
| errorcall(call, _("non-conformable arguments")); |
| } |
| else if (PRIMVAL(op) == 1) { |
| if (nrx != nry) |
| error(_("non-conformable arguments")); |
| } |
| else { |
| if (ncx != ncy) |
| error(_("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, dnx=R_NilValue, dny=R_NilValue; |
| |
| /* allocate dimnames and dimnamesnames */ |
| |
| PROTECT(dimnames = allocVector(VECSXP, 2)); |
| PROTECT(dimnamesnames = allocVector(STRSXP, 2)); |
| if (xdims != R_NilValue) { |
| if (ldx == 2 || ncx == 1) { |
| SET_VECTOR_ELT(dimnames, 0, VECTOR_ELT(xdims, 0)); |
| dnx = getAttrib(xdims, R_NamesSymbol); |
| if(!isNull(dnx)) |
| SET_STRING_ELT(dimnamesnames, 0, STRING_ELT(dnx, 0)); |
| } |
| } |
| |
| #define YDIMS_ET_CETERA \ |
| if (ydims != R_NilValue) { \ |
| if (ldy == 2) { \ |
| SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 1)); \ |
| dny = getAttrib(ydims, R_NamesSymbol); \ |
| if(!isNull(dny)) \ |
| SET_STRING_ELT(dimnamesnames, 1, STRING_ELT(dny, 1)); \ |
| } else if (nry == 1) { \ |
| SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 0)); \ |
| dny = getAttrib(ydims, R_NamesSymbol); \ |
| if(!isNull(dny)) \ |
| SET_STRING_ELT(dimnamesnames, 1, STRING_ELT(dny, 0)); \ |
| } \ |
| } \ |
| \ |
| /* We sometimes attach a dimnames attribute \ |
| * whose elements are all NULL ... \ |
| * This is ugly but causes no real damage. \ |
| * Now (2.1.0 ff), we don't anymore: */ \ |
| if (VECTOR_ELT(dimnames,0) != R_NilValue || \ |
| VECTOR_ELT(dimnames,1) != R_NilValue) { \ |
| if (dnx != R_NilValue || dny != R_NilValue) \ |
| setAttrib(dimnames, R_NamesSymbol, dimnamesnames); \ |
| setAttrib(ans, R_DimNamesSymbol, dimnames); \ |
| } \ |
| UNPROTECT(2) |
| |
| YDIMS_ET_CETERA; |
| } |
| } |
| |
| else if (PRIMVAL(op) == 1) { /* 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)); |
| |
| if (xdims != R_NilValue) { |
| if (ldx == 2) {/* not nrx==1 : .. fixed, ihaka 2003-09-30 */ |
| SET_VECTOR_ELT(dimnames, 0, VECTOR_ELT(xdims, 1)); |
| dnx = getAttrib(xdims, R_NamesSymbol); |
| if(!isNull(dnx)) |
| SET_STRING_ELT(dimnamesnames, 0, STRING_ELT(dnx, 1)); |
| } |
| } |
| |
| YDIMS_ET_CETERA; |
| } |
| |
| } |
| else { /* op == 2: tcrossprod() */ |
| |
| PROTECT(ans = allocMatrix(mode, nrx, nry)); |
| if (mode == CPLXSXP) |
| if(sym) |
| tccrossprod(COMPLEX(CAR(args)), nrx, ncx, |
| COMPLEX(CAR(args)), nry, ncy, COMPLEX(ans)); |
| else |
| tccrossprod(COMPLEX(CAR(args)), nrx, ncx, |
| COMPLEX(CADR(args)), nry, ncy, COMPLEX(ans)); |
| else { |
| if(sym) |
| symtcrossprod(REAL(CAR(args)), nrx, ncx, REAL(ans)); |
| else |
| tcrossprod(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)); |
| |
| if (xdims != R_NilValue) { |
| if (ldx == 2) { |
| SET_VECTOR_ELT(dimnames, 0, VECTOR_ELT(xdims, 0)); |
| dnx = getAttrib(xdims, R_NamesSymbol); |
| if(!isNull(dnx)) |
| SET_STRING_ELT(dimnamesnames, 0, STRING_ELT(dnx, 0)); |
| } |
| } |
| if (ydims != R_NilValue) { |
| if (ldy == 2) { |
| SET_VECTOR_ELT(dimnames, 1, VECTOR_ELT(ydims, 0)); |
| dny = getAttrib(ydims, R_NamesSymbol); |
| if(!isNull(dny)) |
| SET_STRING_ELT(dimnamesnames, 1, STRING_ELT(dny, 0)); |
| } |
| } |
| if (VECTOR_ELT(dimnames,0) != R_NilValue || |
| VECTOR_ELT(dimnames,1) != R_NilValue) { |
| if (dnx != R_NilValue || dny != R_NilValue) |
| setAttrib(dimnames, R_NamesSymbol, dimnamesnames); |
| setAttrib(ans, R_DimNamesSymbol, dimnames); |
| } |
| |
| UNPROTECT(2); |
| } |
| } |
| UNPROTECT(3); |
| return ans; |
| } |
| #undef YDIMS_ET_CETERA |
| |
| SEXP attribute_hidden do_transpose(SEXP call, SEXP op, SEXP args, SEXP rho) |
| { |
| SEXP a, r, dims, dimnames, dimnamesnames = R_NilValue, |
| ndimnamesnames, rnames, cnames; |
| int ldim, ncol = 0, nrow = 0; |
| R_xlen_t len = 0; |
| |
| checkArity(op, args); |
| a = CAR(args); |
| |
| if (isVector(a)) { |
| dims = getAttrib(a, R_DimSymbol); |
| ldim = length(dims); |
| rnames = R_NilValue; |
| cnames = R_NilValue; |
| switch(ldim) { |
| case 0: |
| len = nrow = LENGTH(a); |
| ncol = 1; |
| rnames = getAttrib(a, R_NamesSymbol); |
| dimnames = rnames;/* for isNull() below*/ |
| break; |
| case 1: |
| len = nrow = LENGTH(a); |
| ncol = 1; |
| dimnames = getAttrib(a, R_DimNamesSymbol); |
| if (dimnames != R_NilValue) { |
| rnames = VECTOR_ELT(dimnames, 0); |
| dimnamesnames = getAttrib(dimnames, R_NamesSymbol); |
| } |
| break; |
| case 2: |
| ncol = ncols(a); |
| nrow = nrows(a); |
| len = XLENGTH(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(dimnamesnames); |
| PROTECT(r = allocVector(TYPEOF(a), len)); |
| R_xlen_t i, j, l_1 = len-1; |
| switch (TYPEOF(a)) { |
| case LGLSXP: |
| case INTSXP: |
| // filling in columnwise, "accessing row-wise": |
| for (i = 0, j = 0; i < len; i++, j += nrow) { |
| if (j > l_1) j -= l_1; |
| INTEGER(r)[i] = INTEGER(a)[j]; |
| } |
| break; |
| case REALSXP: |
| for (i = 0, j = 0; i < len; i++, j += nrow) { |
| if (j > l_1) j -= l_1; |
| REAL(r)[i] = REAL(a)[j]; |
| } |
| break; |
| case CPLXSXP: |
| for (i = 0, j = 0; i < len; i++, j += nrow) { |
| if (j > l_1) j -= l_1; |
| COMPLEX(r)[i] = COMPLEX(a)[j]; |
| } |
| break; |
| case STRSXP: |
| for (i = 0, j = 0; i < len; i++, j += nrow) { |
| if (j > l_1) j -= l_1; |
| SET_STRING_ELT(r, i, STRING_ELT(a,j)); |
| } |
| break; |
| case VECSXP: |
| for (i = 0, j = 0; i < len; i++, j += nrow) { |
| if (j > l_1) j -= l_1; |
| SET_VECTOR_ELT(r, i, VECTOR_ELT(a,j)); |
| } |
| break; |
| case RAWSXP: |
| for (i = 0, j = 0; i < len; i++, j += nrow) { |
| if (j > l_1) j -= l_1; |
| RAW(r)[i] = RAW(a)[j]; |
| } |
| break; |
| default: |
| UNPROTECT(2); /* r, dimnamesnames */ |
| goto not_matrix; |
| } |
| PROTECT(dims = allocVector(INTSXP, 2)); |
| INTEGER(dims)[0] = ncol; |
| INTEGER(dims)[1] = nrow; |
| setAttrib(r, R_DimSymbol, dims); |
| UNPROTECT(1); /* dims */ |
| /* R <= 2.2.0: dropped list(NULL,NULL) dimnames : |
| * if(rnames != R_NilValue || cnames != R_NilValue) */ |
| if(!isNull(dimnames)) { |
| 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_VECTOR_ELT(ndimnamesnames, 1, STRING_ELT(dimnamesnames, 0)); |
| SET_VECTOR_ELT(ndimnamesnames, 0, |
| (ldim == 2) ? STRING_ELT(dimnamesnames, 1): |
| R_BlankString); |
| setAttrib(dimnames, R_NamesSymbol, ndimnamesnames); |
| UNPROTECT(1); /* ndimnamesnames */ |
| } |
| setAttrib(r, R_DimNamesSymbol, dimnames); |
| UNPROTECT(1); /* dimnames */ |
| } |
| copyMostAttrib(a, r); |
| UNPROTECT(2); /* r, dimnamesnames */ |
| return r; |
| not_matrix: |
| error(_("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] == isr[itmp]-1) iip[itmp] = 0; \ |
| else { \ |
| iip[itmp]++; \ |
| break; \ |
| } \ |
| for (lj = 0, itmp = 0; itmp < n; itmp++) \ |
| lj += iip[itmp] * stride[itmp]; |
| |
| /* aperm (a, perm, resize = TRUE) */ |
| SEXP attribute_hidden do_aperm(SEXP call, SEXP op, SEXP args, SEXP rho) |
| { |
| SEXP a, perm, r, dimsa, dimsr, dna; |
| int i, j, n, itmp; |
| |
| checkArity(op, args); |
| |
| a = CAR(args); |
| if (!isArray(a)) |
| error(_("invalid first argument, must be an array")); |
| |
| PROTECT(dimsa = getAttrib(a, R_DimSymbol)); |
| n = LENGTH(dimsa); |
| int *isa = INTEGER(dimsa); |
| |
| /* check the permutation */ |
| |
| int *pp = (int *) R_alloc((size_t) n, sizeof(int)); |
| perm = CADR(args); |
| if (length(perm) == 0) { |
| for (i = 0; i < n; i++) pp[i] = n-1-i; |
| } else { |
| if (LENGTH(perm) != n) |
| error(_("'perm' is of wrong length %d (!= %d)"), |
| LENGTH(perm), n); |
| if (isString(perm)) { |
| SEXP dna = getAttrib(a, R_DimNamesSymbol); |
| if (isNull(dna)) |
| error(_("'a' does not have named dimnames")); |
| SEXP dnna = getAttrib(dna, R_NamesSymbol); |
| if (isNull(dnna)) |
| error(_("'a' does not have named dimnames")); |
| for (i = 0; i < n; i++) { |
| const char *this = translateChar(STRING_ELT(perm, i)); |
| for (j = 0; j < n; j++) |
| if (streql(translateChar(STRING_ELT(dnna, j)), |
| this)) {pp[i] = j; break;} |
| if (j >= n) |
| error(_("'perm[%d]' does not match a dimension name"), i+1); |
| } |
| } else { |
| PROTECT(perm = coerceVector(perm, INTSXP)); |
| for (i = 0; i < n; i++) pp[i] = INTEGER(perm)[i] - 1; |
| UNPROTECT(1); |
| } |
| } |
| |
| R_xlen_t *iip = (R_xlen_t *) R_alloc((size_t) n, sizeof(R_xlen_t)); |
| 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 error(_("value out of range in 'perm'")); |
| for (i = 0; i < n; i++) |
| if (iip[i] == 0) error(_("invalid '%s' argument"), "perm"); |
| |
| /* create the stride object and permute */ |
| |
| R_xlen_t *stride = (R_xlen_t *) R_alloc((size_t) n, sizeof(R_xlen_t)); |
| for (iip[0] = 1, i = 1; i<n; i++) iip[i] = iip[i-1] * isa[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)); |
| int *isr = INTEGER(dimsr); |
| for (i = 0; i < n; i++) isr[i] = isa[pp[i]]; |
| |
| /* and away we go! iip will hold the incrementer */ |
| |
| R_xlen_t len = XLENGTH(a); |
| PROTECT(r = allocVector(TYPEOF(a), len)); |
| |
| for (i = 0; i < n; iip[i++] = 0); |
| |
| R_xlen_t li, lj; |
| switch (TYPEOF(a)) { |
| |
| case INTSXP: |
| for (lj = 0, li = 0; li < len; li++) { |
| INTEGER(r)[li] = INTEGER(a)[lj]; |
| CLICKJ; |
| } |
| break; |
| |
| case LGLSXP: |
| for (lj = 0, li = 0; li < len; li++) { |
| LOGICAL(r)[li] = LOGICAL(a)[lj]; |
| CLICKJ; |
| } |
| break; |
| |
| case REALSXP: |
| for (lj = 0, li = 0; li < len; li++) { |
| REAL(r)[li] = REAL(a)[lj]; |
| CLICKJ; |
| } |
| break; |
| |
| case CPLXSXP: |
| for (lj = 0, li = 0; li < len; li++) { |
| COMPLEX(r)[li].r = COMPLEX(a)[lj].r; |
| COMPLEX(r)[li].i = COMPLEX(a)[lj].i; |
| CLICKJ; |
| } |
| break; |
| |
| case STRSXP: |
| for (lj = 0, li = 0; li < len; li++) { |
| SET_STRING_ELT(r, li, STRING_ELT(a, lj)); |
| CLICKJ; |
| } |
| break; |
| |
| case VECSXP: |
| for (lj = 0, li = 0; li < len; li++) { |
| SET_VECTOR_ELT(r, li, VECTOR_ELT(a, lj)); |
| CLICKJ; |
| } |
| break; |
| |
| case RAWSXP: |
| for (lj = 0, li = 0; li < len; li++) { |
| RAW(r)[li] = RAW(a)[lj]; |
| CLICKJ; |
| } |
| break; |
| |
| default: |
| UNIMPLEMENTED_TYPE("aperm", a); |
| } |
| |
| /* handle the resize */ |
| int resize = asLogical(CADDR(args)); |
| if (resize == NA_LOGICAL) error(_("'resize' must be TRUE or FALSE")); |
| |
| /* and handle names(dim(.)) and the dimnames if any */ |
| if (resize) { |
| SEXP nmdm = getAttrib(dimsa, R_NamesSymbol); |
| if(nmdm != R_NilValue) { // dimsr needs correctly permuted names() |
| PROTECT(nmdm); |
| SEXP nm_dr = PROTECT(allocVector(STRSXP, n)); |
| for (i = 0; i < n; i++) { |
| SET_STRING_ELT(nm_dr, i, STRING_ELT(nmdm, pp[i])); |
| } |
| setAttrib(dimsr, R_NamesSymbol, nm_dr); |
| UNPROTECT(2); |
| } |
| setAttrib(r, R_DimSymbol, dimsr); |
| |
| PROTECT(dna = getAttrib(a, R_DimNamesSymbol)); |
| if (dna != R_NilValue) { |
| SEXP dnna, dnr, dnnr; |
| |
| PROTECT(dnr = allocVector(VECSXP, n)); |
| PROTECT(dnna = getAttrib(dna, R_NamesSymbol)); |
| if (dnna != R_NilValue) { |
| PROTECT(dnnr = allocVector(STRSXP, n)); |
| for (i = 0; i < n; i++) { |
| SET_VECTOR_ELT(dnr, i, VECTOR_ELT(dna, pp[i])); |
| SET_STRING_ELT(dnnr, i, STRING_ELT(dnna, pp[i])); |
| } |
| setAttrib(dnr, R_NamesSymbol, dnnr); |
| UNPROTECT(1); |
| } else { |
| for (i = 0; i < n; i++) |
| SET_VECTOR_ELT(dnr, i, VECTOR_ELT(dna, pp[i])); |
| } |
| setAttrib(r, R_DimNamesSymbol, dnr); |
| UNPROTECT(2); |
| } |
| UNPROTECT(1); |
| } |
| else // !resize |
| setAttrib(r, R_DimSymbol, dimsa); |
| |
| UNPROTECT(3); /* dimsa, r, dimsr */ |
| return r; |
| } |
| |
| /* colSums(x, n, p, na.rm) and friends */ |
| SEXP attribute_hidden do_colsum(SEXP call, SEXP op, SEXP args, SEXP rho) |
| { |
| SEXP x, ans = R_NilValue; |
| int type; |
| Rboolean NaRm, keepNA; |
| |
| checkArity(op, args); |
| x = CAR(args); args = CDR(args); |
| R_xlen_t n = asVecSize(CAR(args)); args = CDR(args); |
| R_xlen_t p = asVecSize(CAR(args)); args = CDR(args); |
| NaRm = asLogical(CAR(args)); |
| if (n == NA_INTEGER || n < 0) |
| error(_("invalid '%s' argument"), "n"); |
| if (p == NA_INTEGER || p < 0) |
| error(_("invalid '%s' argument"), "p"); |
| if (NaRm == NA_LOGICAL) error(_("invalid '%s' argument"), "na.rm"); |
| keepNA = !NaRm; |
| |
| switch (type = TYPEOF(x)) { |
| case LGLSXP: |
| case INTSXP: |
| case REALSXP: break; |
| default: |
| error(_("'x' must be numeric")); |
| } |
| if (n * (double)p > XLENGTH(x)) |
| error(_("'x' is too short")); /* PR#16367 */ |
| |
| int OP = PRIMVAL(op); |
| if (OP == 0 || OP == 1) { /* columns */ |
| PROTECT(ans = allocVector(REALSXP, p)); |
| #ifdef _OPENMP |
| int nthreads; |
| /* This gives a spurious -Wunused-but-set-variable error */ |
| if (R_num_math_threads > 0) |
| nthreads = R_num_math_threads; |
| else |
| nthreads = 1; /* for now */ |
| #pragma omp parallel for num_threads(nthreads) default(none) \ |
| firstprivate(x, ans, n, p, type, NaRm, keepNA, R_NaReal, R_NaInt, OP) |
| #endif |
| for (R_xlen_t j = 0; j < p; j++) { |
| R_xlen_t cnt = n, i; |
| LDOUBLE sum = 0.0; |
| switch (type) { |
| case REALSXP: |
| { |
| double *rx = REAL(x) + (R_xlen_t)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;} // unused |
| } |
| break; |
| } |
| case INTSXP: |
| { |
| int *ix = INTEGER(x) + (R_xlen_t)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: |
| { |
| int *ix = LOGICAL(x) + (R_xlen_t)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) sum /= cnt; /* gives NaN for cnt = 0 */ |
| REAL(ans)[j] = (double) sum; |
| } |
| } |
| else { /* rows */ |
| PROTECT(ans = allocVector(REALSXP, n)); |
| |
| /* allocate scratch storage to allow accumulating by columns |
| to improve cache hits */ |
| int *Cnt = NULL; |
| LDOUBLE *rans; |
| if(n <= 10000) { |
| R_CheckStack2(n * sizeof(LDOUBLE)); |
| rans = (LDOUBLE *) alloca(n * sizeof(LDOUBLE)); |
| Memzero(rans, n); |
| } else rans = Calloc(n, LDOUBLE); |
| if (!keepNA && OP == 3) Cnt = Calloc(n, int); |
| |
| for (R_xlen_t j = 0; j < p; j++) { |
| LDOUBLE *ra = rans; |
| switch (type) { |
| case REALSXP: |
| { |
| double *rx = REAL(x) + (R_xlen_t)n * j; |
| if (keepNA) |
| for (R_xlen_t i = 0; i < n; i++) *ra++ += *rx++; |
| else |
| for (R_xlen_t i = 0; i < n; i++, ra++, rx++) |
| if (!ISNAN(*rx)) { |
| *ra += *rx; |
| if (OP == 3) Cnt[i]++; |
| } |
| break; |
| } |
| case INTSXP: |
| { |
| int *ix = INTEGER(x) + (R_xlen_t)n * j; |
| for (R_xlen_t i = 0; i < n; i++, ra++, ix++) |
| if (keepNA) { |
| if (*ix != NA_INTEGER) *ra += *ix; |
| else *ra = NA_REAL; |
| } |
| else if (*ix != NA_INTEGER) { |
| *ra += *ix; |
| if (OP == 3) Cnt[i]++; |
| } |
| break; |
| } |
| case LGLSXP: |
| { |
| int *ix = LOGICAL(x) + (R_xlen_t)n * j; |
| for (R_xlen_t i = 0; i < n; i++, ra++, ix++) |
| if (keepNA) { |
| if (*ix != NA_LOGICAL) *ra += *ix; |
| else *ra = NA_REAL; |
| } |
| else if (*ix != NA_LOGICAL) { |
| *ra += *ix; |
| if (OP == 3) Cnt[i]++; |
| } |
| break; |
| } |
| } |
| } |
| if (OP == 3) { |
| if (keepNA) |
| for (R_xlen_t i = 0; i < n; i++) rans[i] /= p; |
| else |
| for (R_xlen_t i = 0; i < n; i++) rans[i] /= Cnt[i]; |
| } |
| for (R_xlen_t i = 0; i < n; i++) REAL(ans)[i] = (double) rans[i]; |
| |
| if (!keepNA && OP == 3) Free(Cnt); |
| if(n > 10000) Free(rans); |
| } |
| |
| UNPROTECT(1); |
| return ans; |
| } |
| |
| /* |
| { |
| data <- as.vector(data) |
| dim <- as.integer(dim) |
| vl <- prod(dim) |
| if (length(data) != vl) { |
| if (vl > .Machine$integer.max) |
| stop("'dim' specifies too large an array") |
| data <- rep(data, length.out = vl) |
| } |
| if (length(dim)) |
| dim(data) <- dim |
| if (is.list(dimnames) && length(dimnames)) |
| dimnames(data) <- dimnames |
| data |
| } |
| */ |
| |
| /* array(data, dim, dimnames) */ |
| SEXP attribute_hidden do_array(SEXP call, SEXP op, SEXP args, SEXP rho) |
| { |
| SEXP vals, ans, dims, dimnames; |
| R_xlen_t lendat, i, nans; |
| |
| checkArity(op, args); |
| vals = CAR(args); |
| /* at least NULL can get here */ |
| switch(TYPEOF(vals)) { |
| case LGLSXP: |
| case INTSXP: |
| case REALSXP: |
| case CPLXSXP: |
| case STRSXP: |
| case RAWSXP: |
| case EXPRSXP: |
| case VECSXP: |
| break; |
| default: |
| error(_("'data' must be of a vector type, was '%s'"), |
| type2char(TYPEOF(vals))); |
| } |
| lendat = XLENGTH(vals); |
| dims = CADR(args); |
| dimnames = CADDR(args); |
| PROTECT(dims = coerceVector(dims, INTSXP)); |
| int nd = LENGTH(dims); |
| if (nd == 0) error(_("'dims' cannot be of length 0")); |
| double d = 1.0; |
| for (int j = 0; j < nd; j++) d *= INTEGER(dims)[j]; |
| #ifndef LONG_VECTOR_SUPPORT |
| if (d > INT_MAX) error(_("too many elements specified")); |
| #endif |
| nans = (R_xlen_t) d; |
| |
| PROTECT(ans = allocVector(TYPEOF(vals), nans)); |
| switch(TYPEOF(vals)) { |
| case LGLSXP: |
| if (nans && lendat) |
| xcopyLogicalWithRecycle(LOGICAL(ans), LOGICAL(vals), 0, nans, |
| lendat); |
| else |
| for (i = 0; i < nans; i++) LOGICAL(ans)[i] = NA_LOGICAL; |
| break; |
| case INTSXP: |
| if (nans && lendat) |
| xcopyIntegerWithRecycle(INTEGER(ans), INTEGER(vals), 0, nans, |
| lendat); |
| else |
| for (i = 0; i < nans; i++) INTEGER(ans)[i] = NA_INTEGER; |
| break; |
| case REALSXP: |
| if (nans && lendat) |
| xcopyRealWithRecycle(REAL(ans), REAL(vals), 0, nans, lendat); |
| else |
| for (i = 0; i < nans; i++) REAL(ans)[i] = NA_REAL; |
| break; |
| case CPLXSXP: |
| if (nans && lendat) |
| xcopyComplexWithRecycle(COMPLEX(ans), COMPLEX(vals), 0, nans, |
| lendat); |
| else { |
| Rcomplex na_cmplx; |
| na_cmplx.r = NA_REAL; |
| na_cmplx.i = 0; |
| for (i = 0; i < nans; i++) COMPLEX(ans)[i] = na_cmplx; |
| } |
| break; |
| case RAWSXP: |
| if (nans && lendat) |
| xcopyRawWithRecycle(RAW(ans), RAW(vals), 0, nans, lendat); |
| else |
| for (i = 0; i < nans; i++) RAW(ans)[i] = 0; |
| break; |
| case STRSXP: |
| if (nans && lendat) |
| xcopyStringWithRecycle(ans, vals, 0, nans, lendat); |
| else |
| for (i = 0; i < nans; i++) SET_STRING_ELT(ans, i, NA_STRING); |
| break; |
| /* Rest are already initialized */ |
| case VECSXP: |
| case EXPRSXP: |
| #ifdef SWITCH_TO_REFCNT |
| if (nans && lendat) |
| xcopyVectorWithRecycle(ans, vals, 0, nans, lendat); |
| #else |
| if (nans && lendat) { |
| /* Need to guard against possible sharing of values under |
| NAMED. This is not needed with reference |
| coutning. (PR#15919) */ |
| Rboolean needsmark = (lendat < nans || MAYBE_REFERENCED(vals)); |
| for (i = 0; i < nans; i++) { |
| SEXP elt = VECTOR_ELT(vals, i % lendat); |
| if (needsmark || MAYBE_REFERENCED(elt)) |
| MARK_NOT_MUTABLE(elt); |
| SET_VECTOR_ELT(ans, i, elt); |
| } |
| } |
| #endif |
| break; |
| default: |
| // excluded above |
| break; |
| } |
| |
| ans = dimgets(ans, dims); |
| if(!isNull(dimnames) && length(dimnames) > 0) |
| ans = dimnamesgets(ans, dimnames); |
| |
| UNPROTECT(2); |
| return ans; |
| } |
| |
| SEXP attribute_hidden do_diag(SEXP call, SEXP op, SEXP args, SEXP rho) |
| { |
| SEXP ans, x, snr, snc; |
| int nr = 1, nc = 1, nprotect = 1; |
| |
| checkArity(op, args); |
| x = CAR(args); |
| snr = CADR(args); |
| snc = CADDR(args); |
| nr = asInteger(snr); |
| if (nr == NA_INTEGER) |
| error(_("invalid 'nrow' value (too large or NA)")); |
| if (nr < 0) |
| error(_("invalid 'nrow' value (< 0)")); |
| nc = asInteger(snc); |
| if (nc == NA_INTEGER) |
| error(_("invalid 'ncol' value (too large or NA)")); |
| if (nc < 0) |
| error(_("invalid 'ncol' value (< 0)")); |
| int mn = (nr < nc) ? nr : nc; |
| if (mn > 0 && length(x) == 0) |
| error(_("'x' must have positive length")); |
| |
| #ifndef LONG_VECTOR_SUPPORT |
| if ((double)nr * (double)nc > INT_MAX) |
| error(_("too many elements specified")); |
| #endif |
| |
| int nx = LENGTH(x); |
| R_xlen_t NR = nr; |
| |
| #define mk_DIAG(_zero_) \ |
| for (R_xlen_t i = 0; i < NR*nc; i++) ra[i] = _zero_; \ |
| R_xlen_t i, i1; \ |
| MOD_ITERATE1(mn, nx, i, i1, { \ |
| ra[i * (NR+1)] = rx[i1]; \ |
| }); |
| |
| switch(TYPEOF(x)) { |
| |
| case REALSXP: |
| { |
| #define mk_REAL_DIAG \ |
| PROTECT(ans = allocMatrix(REALSXP, nr, nc)); \ |
| double *rx = REAL(x), *ra = REAL(ans); \ |
| mk_DIAG(0.0) |
| |
| mk_REAL_DIAG; |
| break; |
| } |
| case CPLXSXP: |
| { |
| PROTECT(ans = allocMatrix(CPLXSXP, nr, nc)); |
| int nx = LENGTH(x); |
| R_xlen_t NR = nr; |
| Rcomplex *rx = COMPLEX(x), *ra = COMPLEX(ans), zero; |
| zero.r = zero.i = 0.0; |
| mk_DIAG(zero); |
| break; |
| } |
| case INTSXP: |
| { |
| PROTECT(ans = allocMatrix(INTSXP, nr, nc)); |
| int *rx = INTEGER(x), *ra = INTEGER(ans); |
| mk_DIAG(0); |
| break; |
| } |
| case LGLSXP: |
| { |
| PROTECT(ans = allocMatrix(LGLSXP, nr, nc)); |
| int *rx = LOGICAL(x), *ra = LOGICAL(ans); |
| mk_DIAG(0); |
| break; |
| } |
| case RAWSXP: |
| { |
| PROTECT(ans = allocMatrix(RAWSXP, nr, nc)); |
| Rbyte *rx = RAW(x), *ra = RAW(ans); |
| mk_DIAG((Rbyte) 0); |
| break; |
| } |
| default: { |
| PROTECT(x = coerceVector(x, REALSXP)); |
| nprotect++; |
| mk_REAL_DIAG; |
| } |
| } |
| #undef mk_REAL_DIAG |
| #undef mk_DIAG |
| UNPROTECT(nprotect); |
| return ans; |
| } |
| |
| |
| /* backsolve(r, b, k, upper.tri, transpose) */ |
| SEXP attribute_hidden do_backsolve(SEXP call, SEXP op, SEXP args, SEXP rho) |
| { |
| int nprot = 1; |
| checkArity(op, args); |
| |
| SEXP r = CAR(args); args = CDR(args); |
| SEXP b = CAR(args); args = CDR(args); |
| int nrr = nrows(r), nrb = nrows(b), ncb = ncols(b); |
| int k = asInteger(CAR(args)); args = CDR(args); |
| /* k is the number of rows to be used: there must be at least that |
| many rows and cols in the rhs and at least that many rows on |
| the rhs. |
| */ |
| if (k == NA_INTEGER || k <= 0 || k > nrr || k > ncols(r) || k > nrb) |
| error(_("invalid '%s' argument"), "k"); |
| int upper = asLogical(CAR(args)); args = CDR(args); |
| if (upper == NA_INTEGER) error(_("invalid '%s' argument"), "upper.tri"); |
| int trans = asLogical(CAR(args)); |
| if (trans == NA_INTEGER) error(_("invalid '%s' argument"), "transpose"); |
| if (TYPEOF(r) != REALSXP) {PROTECT(r = coerceVector(r, REALSXP)); nprot++;} |
| if (TYPEOF(b) != REALSXP) {PROTECT(b = coerceVector(b, REALSXP)); nprot++;} |
| double *rr = REAL(r); |
| |
| /* check for zeros on diagonal of r: only k row/cols are used. */ |
| size_t incr = nrr + 1; |
| for(int i = 0; i < k; i++) { /* check for zeros on diagonal */ |
| if (rr[i * incr] == 0.0) |
| error(_("singular matrix in 'backsolve'. First zero in diagonal [%d]"), |
| i + 1); |
| } |
| |
| SEXP ans = PROTECT(allocMatrix(REALSXP, k, ncb)); |
| if (k > 0 && ncb > 0) { |
| /* copy (part) cols of b to ans */ |
| for(R_xlen_t j = 0; j < ncb; j++) |
| memcpy(REAL(ans) + j*k, REAL(b) + j*nrb, (size_t)k *sizeof(double)); |
| double one = 1.0; |
| F77_CALL(dtrsm)("L", upper ? "U" : "L", trans ? "T" : "N", "N", |
| &k, &ncb, &one, rr, &nrr, REAL(ans), &k |
| FCONE FCONE FCONE FCONE); |
| } |
| UNPROTECT(nprot); |
| return ans; |
| } |
| |
| /* max.col(m, ties.method) */ |
| SEXP attribute_hidden do_maxcol(SEXP call, SEXP op, SEXP args, SEXP rho) |
| { |
| checkArity(op, args); |
| SEXP m = CAR(args); |
| int method = asInteger(CADR(args)); |
| int nr = nrows(m), nc = ncols(m), nprot = 1; |
| if (TYPEOF(m) != REALSXP) {PROTECT(m = coerceVector(m, REALSXP)); nprot++;} |
| SEXP ans = PROTECT(allocVector(INTSXP, nr)); |
| R_max_col(REAL(m), &nr, &nc, INTEGER(ans), &method); |
| UNPROTECT(nprot); |
| return ans; |
| } |