The R Project SVN R

Rev

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

Rev 19500 Rev 22582
Line 14... Line 14...
14 
#define RELTOL 1e-5
 14 
#define RELTOL 1e-5
15 
 
 15 
 
16 
void R_max_col(double *matrix, int *nr, int *nc, int *maxes)
 16 
void R_max_col(double *matrix, int *nr, int *nc, int *maxes)
17 
{
 17 
{
18 
    int	  r, c, m, ntie, n_r = *nr;
 18 
    int	  r, c, m, ntie, n_r = *nr;
19 
    double a, b, tol;
 19 
    double a, b, tol, large;
- 
 
 20 
    Rboolean isna, used_random=FALSE;
20 
 
 21 
 
- 
 
 22 
    for (r = 0; r < n_r; r++) {
- 
 
 23 
	/* first check row for any NAs and find the largest entry */
- 
 
 24 
	large = 0.0;
21 
#define NA_check				\
 25 
	isna = FALSE;
22 
	if (ISNAN(a)) { /* NA or NaN */		\
 26 
	for (c = 0; c < *nc; c++) {
23 
	    maxes[r] = NA_INTEGER;		\
 27 
	    a = matrix[r + c * n_r];
- 
 
 28 
	    if (ISNAN(a)) { isna = TRUE; break; }
24 
	    continue;				\
 29 
	    large = fmax2(large, fabs(a));
25 
	}
 30 
	}
- 
 
 31 
	if (isna) { maxes[r] = NA_INTEGER; continue; }
- 
 
 32 
	tol = RELTOL * large;
26 
 
 33 
 
27 
    GetRNGstate();
 - 
 
28 
    for (r = 0; r < n_r; r++) {
 - 
 
29 
	m = 0;
 34 
	m = 0;
30 
	ntie = 1;
 35 
	ntie = 1;
31 
	a = matrix[r];
 36 
	a = matrix[r];
32 
	NA_check;
 - 
 
33 
	for (c = 1; c < *nc; c++) {
 37 
	for (c = 1; c < *nc; c++) {
34 
	    b = matrix[r + c * n_r];
 38 
	    b = matrix[r + c * n_r];
35 
	    NA_check;
 - 
 
36 
	    tol = RELTOL * fmax2(fabs(a), fabs(b));
 - 
 
37 
	    if (b >= a + tol) {
 39 
	    if (b >= a + tol) {
38 
		ntie = 1;
 40 
		ntie = 1;
39 
		a = b;
 41 
		a = b;
40 
		m = c;
 42 
		m = c;
41 
	    } else if (b >= a - tol) {
 43 
	    } else if (b >= a - tol) {
42 
		ntie++;
 44 
		ntie++;
- 
 
 45 
		if (!used_random) { GetRNGstate(); used_random = TRUE; }
43 
		if (ntie * unif_rand() < 1.) m = c;
 46 
		if (ntie * unif_rand() < 1.) m = c;
44 
	    }
 47 
	    }
45 
	}
 48 
	}
46 
	maxes[r] = m + 1;
 49 
	maxes[r] = m + 1;
47 
    }
 50 
    }
48 
    PutRNGstate();
 51 
    if(used_random) PutRNGstate();
49 
}
 52 
}