The R Project SVN R

Rev

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

Rev 35275 Rev 35296
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>     /* for Calloc/Free */
31 
#include <R_ext/Applic.h> /* for dgemm */
 32 
#include <R_ext/Applic.h> /* for dgemm */
32 
 
 33 
 
33 
/* "GetRowNames" and "GetColNames" are utility routines which
 34 
/* "GetRowNames" and "GetColNames" are utility routines which
34 
 * locate and return the row names and column names from the
 35 
 * locate and return the row names and column names from the
35 
 * dimnames attribute of a matrix.  They are useful because
 36 
 * dimnames attribute of a matrix.  They are useful because
Line 1084... Line 1085...
1084 
	/* reverse summation order to improve cache hits */
 1085 
	/* reverse summation order to improve cache hits */
1085 
	if (type == REALSXP) {
 1086 
	if (type == REALSXP) {
1086 
	    double *rans = REAL(ans), *ra = rans;
 1087 
	    double *rans = REAL(ans), *ra = rans;
1087 
	    int *Cnt = NULL, *c;
 1088 
	    int *Cnt = NULL, *c;
1088 
	    rx = REAL(x);
 1089 
	    rx = REAL(x);
1089 
	    if (!keepNA && OP == 3) {
 1090 
	    if (!keepNA && OP == 3) Cnt = Calloc(n, int);
1090 
		Cnt = (int *) alloca(n*sizeof(int));
 - 
 
1091 
		memset(Cnt, 0, n*sizeof(int));
 - 
 
1092 
	    }
 - 
 
1093 
	    memset(rans, 0, n*sizeof(double));
 1091 
	    memset(rans, 0, n*sizeof(double));
1094 
	    for (j = 0; j < p; j++) {
 1092 
	    for (j = 0; j < p; j++) {
1095 
		ra = rans;
 1093 
		ra = rans;
1096 
		if (keepNA)
 1094 
		if (keepNA)
1097 
		    for (i = 0; i < n; i++) *ra++ += *rx++;
 1095 
		    for (i = 0; i < n; i++) *ra++ += *rx++;
Line 1102... Line 1100...
1102 
			    if (OP == 3) (*c)++;
 1100 
			    if (OP == 3) (*c)++;
1103 
			}
 1101 
			}
1104 
	    }
 1102 
	    }
1105 
	    if (OP == 3) {
 1103 
	    if (OP == 3) {
1106 
		if (keepNA)
 1104 
		if (keepNA)
1107 
		    for (ra = rans, i = 0; i < n; i++)
 1105 
		    for (ra = rans, i = 0; i < n; i++) *ra++ /= p;
1108 
			*ra++ /= p;
 - 
 
1109 
		else {
 1106 
		else {
1110 
		    for (ra = rans, c = Cnt, i = 0; i < n; i++, c++)
 1107 
		    for (ra = rans, c = Cnt, i = 0; i < n; i++, c++)
1111 
			if (*c > 0) *ra++ /= *c; else *ra++ = NA_REAL;
 1108 
			if (*c > 0) *ra++ /= *c; else *ra++ = NA_REAL;
- 
 
 1109 
		    Free(Cnt);
1112 
		}
 1110 
		}
1113 
	    }
 1111 
	    }
1114 
	    UNPROTECT(1);
 1112 
	    UNPROTECT(1);
1115 
	    return ans;
 1113 
	    return ans;
1116 
	}
 1114 
	}