Google Git
Sign in
third-party-mirror / R / refs/heads/R-3-6-3-tag / . / src / main / array.c
blob: 182aa18f44899c05d98aeb3452336196715afde8 [file] [log] [blame]
/*
* 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;
}