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



#include <config.h>
#endif
 
#include <Defn.h>
#include <Rmath.h>
#include <R_ext/RS.h>     /* for Calloc/Free */
#include <R_ext/Applic.h> /* for dgemm */
 
/* "GetRowNames" and "GetColNames" are utility routines which
 * locate and return the row names and column names from the
 * dimnames attribute of a matrix.  They are useful because

			    x, &nrx, y, &nry, &zero, z, &nrx);
    } else /* zero-extent operations should return zeroes */
	for(i = 0; i < nrx*ncy; i++) z[i] = 0;
}
 
#ifdef HAVE_FORTRAN_DOUBLE_COMPLEX
/* ZGEMM - perform one of the matrix-matrix operations    */
/* C := alpha*op( A )*op( B ) + beta*C */
extern void
F77_NAME(zgemm)(const char *transa, const char *transb, const int *m,
		const int *n, const int *k, const Rcomplex *alpha,

#endif
 
static void cmatprod(Rcomplex *x, int nrx, int ncx,
		     Rcomplex *y, int nry, int ncy, Rcomplex *z)
{
#ifdef HAVE_FORTRAN_DOUBLE_COMPLEX
    char *transa = "N", *transb = "N";
    int i;
    Rcomplex one, zero;
 
    one.r = 1.0; one.i = zero.r = zero.i = 0.0;

	cnt = p;
	PROTECT(ans = allocVector(REALSXP, n));
 
	/* reverse summation order to improve cache hits */
	if (type == REALSXP) {
	    double *rans = REAL(ans), *ra = rans;
	    int *Cnt = NULL, *c;
	    rx = REAL(x);
	    if (!keepNA && OP == 3) Cnt = Calloc(n, int);
	    memset(rans, 0, n*sizeof(double));
	    for (j = 0; j < p; j++) {
		ra = rans;
		if (keepNA)
		    for (i = 0; i < n; i++) *ra++ += *rx++;
		else

			    if (OP == 3) (*c)++;
			}
	    }
	    if (OP == 3) {
		if (keepNA)
		    for (ra = rans, i = 0; i < n; i++) *ra++ /= p;
 
		else {
		    for (ra = rans, c = Cnt, i = 0; i < n; i++, c++)
			if (*c > 0) *ra++ /= *c; else *ra++ = NA_REAL;
		    Free(Cnt);
		}

Rev 34864	Rev 35372
Line 26...	Line 26...
26	`#include <config.h>`	26	`#include <config.h>`
27	`#endif`	27	`#endif`
28		28
29	`#include <Defn.h>`	29	`#include <Defn.h>`
30	`#include <Rmath.h>`	30	`#include <Rmath.h>`
31	`#include <R_ext/RS.h>`	31	`#include <R_ext/RS.h> /* for Calloc/Free */`
32	`#include <R_ext/Applic.h> /* for dgemm */`	32	`#include <R_ext/Applic.h> /* for dgemm */`
33		33
34	`/* "GetRowNames" and "GetColNames" are utility routines which`	34	`/* "GetRowNames" and "GetColNames" are utility routines which`
35	`* locate and return the row names and column names from the`	35	`* locate and return the row names and column names from the`
36	`* dimnames attribute of a matrix. They are useful because`	36	`* dimnames attribute of a matrix. They are useful because`
Line 413...	Line 413...
413	`x, &nrx, y, &nry, &zero, z, &nrx);`	413	`x, &nrx, y, &nry, &zero, z, &nrx);`
414	`} else /* zero-extent operations should return zeroes */`	414	`} else /* zero-extent operations should return zeroes */`
415	`for(i = 0; i < nrx*ncy; i++) z[i] = 0;`	415	`for(i = 0; i < nrx*ncy; i++) z[i] = 0;`
416	`}`	416	`}`
417		417
418	`#ifdef HAVE_DOUBLE_COMPLEX`	418	`#ifdef HAVE_FORTRAN_DOUBLE_COMPLEX`
419	`/* ZGEMM - perform one of the matrix-matrix operations */`	419	`/* ZGEMM - perform one of the matrix-matrix operations */`
420	`/* C := alphaop( A )op( B ) + betaC /`	420	`/* C := alphaop( A )op( B ) + betaC /`
421	`extern void`	421	`extern void`
422	`F77_NAME(zgemm)(const char transa, const char transb, const int *m,`	422	`F77_NAME(zgemm)(const char transa, const char transb, const int *m,`
423	`const int n, const int k, const Rcomplex *alpha,`	423	`const int n, const int k, const Rcomplex *alpha,`
Line 427...	Line 427...
427	`#endif`	427	`#endif`
428		428
429	`static void cmatprod(Rcomplex *x, int nrx, int ncx,`	429	`static void cmatprod(Rcomplex *x, int nrx, int ncx,`
430	`Rcomplex y, int nry, int ncy, Rcomplex z)`	430	`Rcomplex y, int nry, int ncy, Rcomplex z)`
431	`{`	431	`{`
432	`#ifdef HAVE_DOUBLE_COMPLEX`	432	`#ifdef HAVE_FORTRAN_DOUBLE_COMPLEX`
433	`char transa = "N", transb = "N";`	433	`char transa = "N", transb = "N";`
434	`int i;`	434	`int i;`
435	`Rcomplex one, zero;`	435	`Rcomplex one, zero;`
436		436
437	`one.r = 1.0; one.i = zero.r = zero.i = 0.0;`	437	`one.r = 1.0; one.i = zero.r = zero.i = 0.0;`
Line 1082...	Line 1082...
1082	`cnt = p;`	1082	`cnt = p;`
1083	`PROTECT(ans = allocVector(REALSXP, n));`	1083	`PROTECT(ans = allocVector(REALSXP, n));`
1084		1084
1085	`/* reverse summation order to improve cache hits */`	1085	`/* reverse summation order to improve cache hits */`
1086	`if (type == REALSXP) {`	1086	`if (type == REALSXP) {`
1087	`double rans = REAL(ans), ra = rans, Cnt = NULL, c;`	1087	`double rans = REAL(ans), ra = rans;`
-		1088	`int Cnt = NULL, c;`
1088	`rx = REAL(x);`	1089	`rx = REAL(x);`
1089	`if (!keepNA && OP == 3) Cnt = Calloc(n, double);`	1090	`if (!keepNA && OP == 3) Cnt = Calloc(n, int);`
1090	`for (ra = rans, i = 0; i < n; i++) *ra++ = 0.0;`	1091	`memset(rans, 0, n*sizeof(double));`
1091	`for (j = 0; j < p; j++) {`	1092	`for (j = 0; j < p; j++) {`
1092	`ra = rans;`	1093	`ra = rans;`
1093	`if (keepNA)`	1094	`if (keepNA)`
1094	`for (i = 0; i < n; i++) ra++ += rx++;`	1095	`for (i = 0; i < n; i++) ra++ += rx++;`
1095	`else`	1096	`else`
Line 1099...	Line 1100...
1099	`if (OP == 3) (*c)++;`	1100	`if (OP == 3) (*c)++;`
1100	`}`	1101	`}`
1101	`}`	1102	`}`
1102	`if (OP == 3) {`	1103	`if (OP == 3) {`
1103	`if (keepNA)`	1104	`if (keepNA)`
1104	`for (ra = rans, i = 0; i < n; i++)`	1105	`for (ra = rans, i = 0; i < n; i++) *ra++ /= p;`
1105	`*ra++ /= p;`	-
1106	`else {`	1106	`else {`
1107	`for (ra = rans, c = Cnt, i = 0; i < n; i++, c++)`	1107	`for (ra = rans, c = Cnt, i = 0; i < n; i++, c++)`
1108	`if (c > 0) ra++ /= c; else ra++ = NA_REAL;`	1108	`if (c > 0) ra++ /= c; else ra++ = NA_REAL;`
1109	`Free(Cnt);`	1109	`Free(Cnt);`
1110	`}`	1110	`}`

The R Project SVN R

(root)/branches/R-noNULLenv/src/main/array.c @ 34412 – Rev 34864 → 35372