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

 #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.
     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 = getElementLength(x, i, call, rho);
+        *ans_elt = (double) getElementLength(x, i, call, rho);
     UNPROTECT(1);
     return ans;
+}
 #endif
 	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 (!R_FINITE(x[i].r+x[i].i+x[i+1].r+x[i+1].i))
 	    return TRUE;
     return FALSE;
+}
-static R_INLINE void simple_matprod(double *x, int nrx, int ncx,
+/* experimental version for SIMD hardware (see also mayHaveNaNOrInf_simd) */
-                                    double *y, int nry, int ncy, double *z)
+static Rboolean cmayHaveNaNOrInf_simd(Rcomplex *x, R_xlen_t n)
+{
-    LDOUBLE sum;
+    double s = 0;
+    /* _OPENMP >= 201307 - see mayHaveNaNOrInf_simd */
+#if defined(_OPENMP) && HAVE_OPENMP_SIMDRED
+    #pragma omp simd reduction(+:s)
+#endif
-    R_xlen_t NRX = nrx, NRY = nry;
+    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 i = 0; i < nrx; i++)			\
-	for (int k = 0; k < ncy; k++) {
+	for (int k = 0; k < ncy; k++) {			\
-	    sum = 0.0;
+	    sum = 0.0;					\
-	    for (int j = 0; j < ncx; j++)
+	    for (int j = 0; j < ncx; j++)		\
-		sum += x[i + j * NRX] * y[j + k * NRY];
+		sum += x[i + j * NRX] * y[j + k * NRY];	\
-	    z[i + k * NRX] = (double) sum;
+	    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 R_INLINE void simple_crossprod(double *x, int nrx, int ncx,
+static void internal_crossprod(double *x, int nrx, int ncx,
-                                      double *y, int nry, int ncy, double *z)
+                               double *y, int nry, int ncy, double *z)
+{
     LDOUBLE sum;
+#define CROSSPROD_BODY					\
-    R_xlen_t NRX = nrx, NRY = nry, NCX = ncx;
+    R_xlen_t NRX = nrx, NRY = nry, NCX = ncx;		\
-    for (int i = 0; i < ncx; i++)
+    for (int i = 0; i < ncx; i++)			\
-	for (int k = 0; k < ncy; k++) {
+	for (int k = 0; k < ncy; k++) {			\
-	    sum = 0.0;
+	    sum = 0.0;					\
-	    for (int j = 0; j < nrx; j++)
+	    for (int j = 0; j < nrx; j++)		\
-		sum += x[j + i * NRX] * y[j + k * NRY];
+		sum += x[j + i * NRX] * y[j + k * NRY];	\
-	    z[i + k * NCX] = (double) sum;
+	    z[i + k * NCX] = (double) sum;		\
+	}
+    CROSSPROD_BODY;
+}
-static R_INLINE void simple_tcrossprod(double *x, int nrx, int ncx,
+static void simple_crossprod(double *x, int nrx, int ncx,
-                                       double *y, int nry, int ncy, double *z)
+                             double *y, int nry, int ncy, double *z)
+{
-    LDOUBLE sum;
+    double sum;
-    R_xlen_t NRX = nrx, NRY = nry;
+    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 i = 0; i < nrx; i++)			\
-	for (int k = 0; k < nry; k++) {
+	for (int k = 0; k < nry; k++) {			\
-	    sum = 0.0;
+	    sum = 0.0;					\
-	    for (int j = 0; j < ncx; j++)
+	    for (int j = 0; j < ncx; j++)		\
-		sum += x[i + j * NRX] * y[k + j * NRY];
+		sum += x[i + j * NRX] * y[k + j * NRY];	\
-	    z[i + k * NRX] = (double) sum;
+	    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)
+{
-    char *transN = "N", *transT = "T";
-    double one = 1.0, zero = 0.0;
-    int ione = 1;
     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;
+    }
-    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
+    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))
+	    if (mayHaveNaNOrInf(x, NRX*ncx) || mayHaveNaNOrInf(y, NRY*ncy)) {
-	    simple_matprod(x, nrx, ncx, y, nry, ncy, z);
+		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;
-	else if (ncy == 1) /* matrix-vector or dot product */
+    if (ncy == 1) /* matrix-vector or dot product */
-	    F77_CALL(dgemv)(transN, &nrx, &ncx, &one, x,
+	F77_CALL(dgemv)(transN, &nrx, &ncx, &one, x,
-	                    &nrx, y, &ione, &zero, z, &ione);
+			&nrx, y, &ione, &zero, z, &ione);
-	else if (nrx == 1) /* vector-matrix */
+    else if (nrx == 1) /* vector-matrix */
-	    /* Instead of xY, compute (xY)^T == (Y^T)(x^T)
+	/* Instead of xY, compute (xY)^T == (Y^T)(x^T)
-	       The result is a vector, so transposing its content is no-op */
+	   The result is a vector, so transposing its content is no-op */
-	    F77_CALL(dgemv)(transT, &nry, &ncy, &one, y,
+	F77_CALL(dgemv)(transT, &nry, &ncy, &one, y,
-	                    &nry, x, &ione, &zero, z, &ione);
+			&nry, x, &ione, &zero, z, &ione);
-	else /* matrix-matrix or outer product */
+    else /* matrix-matrix or outer product */
-	    F77_CALL(dgemm)(transN, transN, &nrx, &ncy, &ncx, &one,
+	F77_CALL(dgemm)(transN, transN, &nrx, &ncy, &ncx, &one,
-			    x, &nrx, y, &nry, &zero, z, &nrx);
+			x, &nrx, y, &nry, &zero, z, &nrx);
-    } else /* zero-extent operations should return zeroes */
-	for(R_xlen_t i = 0; i < NRX*ncy; i++) z[i] = 0;
+}
-static R_INLINE
-void simple_cmatprod(Rcomplex *x, int nrx, int ncx,
+static void internal_cmatprod(Rcomplex *x, int nrx, int ncx,
-		     Rcomplex *y, int nry, int ncy, Rcomplex *z)
+                              Rcomplex *y, int nry, int ncy, Rcomplex *z)
+{
-    int i, j, k;
-    double xij_r, xij_i, yjk_r, yjk_i;
     LDOUBLE sum_i, sum_r;
+#define CMATPROD_BODY					    \
+    int i, j, k;					    \
+    double complex xij, yjk;				    \
-    R_xlen_t NRX = nrx, NRY = nry;
+    R_xlen_t NRX = nrx, NRY = nry;			    \
-    for (i = 0; i < nrx; i++)
+    for (i = 0; i < nrx; i++)				    \
-	for (k = 0; k < ncy; k++) {
+	for (k = 0; k < ncy; k++) {			    \
-	    sum_r = 0.0;
+	    sum_r = 0.0;				    \
-	    sum_i = 0.0;
+	    sum_i = 0.0;				    \
-	    for (j = 0; j < ncx; j++) {
+	    for (j = 0; j < ncx; j++) {			    \
-		xij_r = x[i + j * NRX].r;
-		xij_i = x[i + j * NRX].i;
+		xij = toC99(x + (i + j * NRX));		    \
-		yjk_r = y[j + k * NRY].r;
-		yjk_i = y[j + k * NRY].i;
+		yjk = toC99(y + (j + k * NRY));		    \
-		sum_r += (xij_r * yjk_r - xij_i * yjk_i);
+		sum_r += creal(xij * yjk);		    \
-		sum_i += (xij_r * yjk_i + xij_i * yjk_r);
+		sum_i += cimag(xij * yjk);		    \
-	    }
+	    }						    \
-	    z[i + k * NRX].r = (double) sum_r;
+	    z[i + k * NRX].r = (double) sum_r;		    \
-	    z[i + k * NRX].i = (double) sum_i;
+	    z[i + k * NRX].i = (double) sum_i;		    \
+	}
+    CMATPROD_BODY;
+}
-static R_INLINE
-void simple_ccrossprod(Rcomplex *x, int nrx, int ncx,
+static void simple_cmatprod(Rcomplex *x, int nrx, int ncx,
-		       Rcomplex *y, int nry, int ncy, Rcomplex *z)
+                            Rcomplex *y, int nry, int ncy, Rcomplex *z)
+{
-    int i, j, k;
-    double xji_r, xji_i, yjk_r, yjk_i;
+    double sum_i, sum_r;
-    LDOUBLE 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;
+    R_xlen_t NRX = nrx, NRY = nry, NCX = ncx;		    \
-    for (i = 0; i < ncx; i++)
+    for (i = 0; i < ncx; i++)				    \
-	for (k = 0; k < ncy; k++) {
+	for (k = 0; k < ncy; k++) {			    \
-	    sum_r = 0.0;
+	    sum_r = 0.0;				    \
-	    sum_i = 0.0;
+	    sum_i = 0.0;				    \
-	    for (j = 0; j < nrx; j++) {
+	    for (j = 0; j < nrx; j++) {			    \
-		xji_r = x[j + i * NRX].r;
-		xji_i = x[j + i * NRX].i;
+		xji = toC99(x + (j + i * NRX));		    \
-		yjk_r = y[j + k * NRY].r;
-		yjk_i = y[j + k * NRY].i;
+		yjk = toC99(y + (j + k * NRY));		    \
-		sum_r += (xji_r * yjk_r - xji_i * yjk_i);
+		sum_r += creal(xji * yjk);		    \
-		sum_i += (xji_r * yjk_i + xji_i * yjk_r);
+		sum_i += cimag(xji * yjk);		    \
-	    }
+	    }						    \
-	    z[i + k * NCX].r = (double) sum_r;
+	    z[i + k * NCX].r = (double) sum_r;		    \
-	    z[i + k * NCX].i = (double) sum_i;
+	    z[i + k * NCX].i = (double) sum_i;		    \
+	}
+    CCROSSPROD_BODY;
+}
-static R_INLINE
-void simple_tccrossprod(Rcomplex *x, int nrx, int ncx,
+static void simple_ccrossprod(Rcomplex *x, int nrx, int ncx,
-		        Rcomplex *y, int nry, int ncy, Rcomplex *z)
+                              Rcomplex *y, int nry, int ncy, Rcomplex *z)
+{
-    int i, j, k;
-    double xij_r, xij_i, ykj_r, ykj_i;
+    double sum_i, sum_r;
-    LDOUBLE 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;
+    R_xlen_t NRX = nrx, NRY = nry;			    \
-    for (i = 0; i < nrx; i++)
+    for (i = 0; i < nrx; i++)				    \
-	for (k = 0; k < nry; k++) {
+	for (k = 0; k < nry; k++) {			    \
-	    sum_r = 0.0;
+	    sum_r = 0.0;				    \
-	    sum_i = 0.0;
+	    sum_i = 0.0;				    \
-	    for (j = 0; j < ncx; j++) {
+	    for (j = 0; j < ncx; j++) {			    \
-		xij_r = x[i + j * NRX].r;
-		xij_i = x[i + j * NRX].i;
+		xij = toC99(x + (i + j * NRX));		    \
-		ykj_r = y[k + j * NRY].r;
-		ykj_i = y[k + j * NRY].i;
+		ykj = toC99(y + (k + j * NRY));		    \
-		sum_r += (xij_r * ykj_r - xij_i * ykj_i);
+		sum_r += creal(xij * ykj);		    \
-		sum_i += (xij_r * ykj_i + xij_i * ykj_r);
+		sum_i += cimag(xij * ykj);		    \
-	    }
+	    }						    \
-	    z[i + k * NRX].r = (double) sum_r;
+	    z[i + k * NRX].r = (double) sum_r;		    \
-	    z[i + k * NRX].i = (double) sum_i;
+	    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)
+{
-    char *transa = "N", *transb = "N";
-    Rcomplex one, zero;
     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;
+    }
-    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
-    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
 #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))
+	    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);
+		simple_cmatprod(x, nrx, ncx, y, nry, ncy, z);
+		return;
-	else
+	    }
+	    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,
+    F77_CALL(zgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
-	                    x, &nrx, y, &nry, &zero, z, &nrx);
+                    x, &nrx, y, &nry, &zero, z, &nrx);
 #endif
-    } else /* zero-extent operations should return zeroes */
-	for(R_xlen_t i = 0; i < NRX*ncy; i++) z[i].r = z[i].i = 0;
+}
 static void symcrossprod(double *x, int nr, int nc, double *z)
+{
-    char *trans = "T", *uplo = "U";
-    double one = 1.0, zero = 0.0;
-    R_xlen_t NC = nc;
+    R_xlen_t NR = nr, NC = nc;
-    if (nr > 0 && nc > 0) {
+    if (nr == 0 || nc == 0) {
-	F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc);
-	for (int i = 1; i < nc; i++)
-	    for (int j = 0; j < i; j++) z[i + NC *j] = z[j + NC * i];
-    } else { /* zero-extent operations should return zeroes */
+	/* 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);
+    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)
+{
-    char *transT = "T", *transN = "N";
-    double one = 1.0, zero = 0.0;
-    int ione = 1;
     R_xlen_t NRX = nrx, NRY = nry;
-    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
+    if (nrx == 0 || ncx == 0 || nry == 0 || ncy == 0) {
-	if (mayHaveNaNOrInf(x, NRX*ncx) || mayHaveNaNOrInf(y, NRY*ncy))
-	    /* see matprod for more details */
-	    simple_crossprod(x, nrx, ncx, y, nry, ncy, z);
-	else if (ncy == 1) /* matrix-vector or dot product */
-	    F77_CALL(dgemv)(transT, &nrx, &ncx, &one, x,
-	                    &nrx, y, &ione, &zero, z, &ione);
-	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);
-	else /* matrix-matrix  or outer product */
-	    F77_CALL(dgemm)(transT, transN, &ncx, &ncy, &nrx, &one,
-			x, &nrx, y, &nry, &zero, z, &ncx);
-    } else { /* zero-extent operations should return zeroes */
+	/* 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);
+    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);
+    else /* matrix-matrix  or outer product */
+	F77_CALL(dgemm)(transT, transN, &ncx, &ncy, &nrx, &one,
+		        x, &nrx, y, &nry, &zero, z, &ncx);
+}
 static void ccrossprod(Rcomplex *x, int nrx, int ncx,
 		       Rcomplex *y, int nry, int ncy, Rcomplex *z)
+{
-    char *transa = "T", *transb = "N";
-    Rcomplex one, zero;
     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;
+    }
-    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
-    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
 #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))
+	    if (cmayHaveNaNOrInf(x, NRX*ncx) || cmayHaveNaNOrInf(y, NRY*ncy)) {
-	    simple_ccrossprod(x, nrx, ncx, y, nry, ncy, z);
+		simple_ccrossprod(x, nrx, ncx, y, nry, ncy, z);
+		return;
-	else
+	    }
-	    F77_CALL(zgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
+	    break; /* use blas */
+	case MATPROD_INTERNAL:
-	                    x, &nrx, y, &nry, &zero, z, &ncx);
+	    internal_ccrossprod(x, nrx, ncx, y, nry, ncy, z);
-#endif
+	    return;
+	case MATPROD_BLAS:
+	    break;
+	case MATPROD_DEFAULT_SIMD:
-    } else { /* zero-extent operations should return zeroes */
+	    if (cmayHaveNaNOrInf_simd(x, NRX*ncx) ||
-	R_xlen_t NCX = ncx;
+		    cmayHaveNaNOrInf_simd(y, NRY*ncy)) {
-	for(R_xlen_t i = 0; i < NCX*ncy; i++) z[i].r = z[i].i = 0;
+		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);
+#endif
+}
 static void symtcrossprod(double *x, int nr, int nc, double *z)
+{
-    char *trans = "N", *uplo = "U";
-    double one = 1.0, zero = 0.0;
+    R_xlen_t NR = nr;
-    if (nr > 0 && nc > 0) {
+    if (nr == 0 || nc == 0) {
-	F77_CALL(dsyrk)(uplo, trans, &nr, &nc, &one, x, &nr, &zero, z, &nr);
-	for (int i = 1; i < nr; i++)
-	    for (int j = 0; j < i; j++) z[i + nr *j] = z[j + nr * i];
-    } else { /* zero-extent operations should return zeroes */
+	/* zero-extent operations should return zeroes */
-	R_xlen_t NR = nr;
 	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);
+    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)
+{
-    char *transN = "N", *transT = "T";
-    double one = 1.0, zero = 0.0;
-    int ione = 1;
     R_xlen_t NRX = nrx, NRY = nry;
-    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
+    if (nrx == 0 || ncx == 0 || nry == 0 || ncy == 0) {
-	if (mayHaveNaNOrInf(x, NRX*ncx) || mayHaveNaNOrInf(y, NRY*ncy))
-	    /* see matprod for more details */
-	    simple_tcrossprod(x, nrx, ncx, y, nry, ncy, z);
-	else if (nry == 1) /* matrix-vector or dot product */
-	    F77_CALL(dgemv)(transN, &nrx, &ncx, &one, x,
-	                    &nrx, y, &ione, &zero, z, &ione);
-	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);
-	else /* matrix-matrix or outer product */
-	    F77_CALL(dgemm)(transN, transT, &nrx, &nry, &ncx, &one,
-			x, &nrx, y, &nry, &zero, z, &nrx);
-    } else { /* zero-extent operations should return zeroes */
+	/* zero-extent operations should return zeroes */
-	R_xlen_t NRX = nrx;
 	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);
+    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);
+    else /* matrix-matrix or outer product */
+	F77_CALL(dgemm)(transN, transT, &nrx, &nry, &ncx, &one,
+		    x, &nrx, y, &nry, &zero, z, &nrx);
+}
 static void tccrossprod(Rcomplex *x, int nrx, int ncx,
 			Rcomplex *y, int nry, int ncy, Rcomplex *z)
+{
-    char *transa = "N", *transb = "T";
-    Rcomplex one, zero;
     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;
+    }
-    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
-    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {
 #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))
+	    if (cmayHaveNaNOrInf(x, NRX*ncx) || cmayHaveNaNOrInf(y, NRY*ncy)) {
-	    simple_tccrossprod(x, nrx, ncx, y, nry, ncy, z);
+		simple_tccrossprod(x, nrx, ncx, y, nry, ncy, z);
+		return;
-	else
+	    }
-	    F77_CALL(zgemm)(transa, transb, &nrx, &nry, &ncx, &one,
+	    break; /* use blas */
+	case MATPROD_INTERNAL:
-	                    x, &nrx, y, &nry, &zero, z, &nrx);
+	    internal_tccrossprod(x, nrx, ncx, y, nry, ncy, z);
-#endif
+	    return;
+	case MATPROD_BLAS:
+	    break;
+	case MATPROD_DEFAULT_SIMD:
-    } else { /* zero-extent operations should return zeroes */
+	    if (cmayHaveNaNOrInf_simd(x, NRX*ncx) ||
-	R_xlen_t NRX = nrx;
+		    cmayHaveNaNOrInf_simd(y, NRY*ncy)) {
-	for(R_xlen_t i = 0; i < NRX*nry; i++) z[i].r = z[i].i = 0;
+		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);
+#endif
+}
 /* "%*%" (op = 0), crossprod (op = 1) or tcrossprod (op = 2) */
 SEXP attribute_hidden do_matprod(SEXP call, SEXP op, SEXP args, SEXP rho)

The R Project SVN R

(root)/branches/ALTREP/src/main/array.c @ 77458 – Rev 72114 → 72834