t_gCMatrix_colSums.c
/*------ Definition of a template for [diln]gCMatrix_colsums(...) : *
* -------- ~~~~~~~~~~~~~~~~~~~~~~
* i.e., included several times from ./dgCMatrix.c
* ~~~~~~~~~~~~~
*/
/* for all cases with an 'x' slot -- i.e. almost all cases ;
* just redefine this in the other cases:
*/
#ifdef _dgC_
# define gCMatrix_colSums dgCMatrix_colSums
# define _DOUBLE_ans
# define _has_x_slot_
/*Future? # define _has_x_d_slot_ */
# undef _dgC_
#elif defined (_igC_)
# define gCMatrix_colSums igCMatrix_colSums
# define _DOUBLE_ans
# define _has_x_slot_
/*Future? # define _has_x_d_slot_ */
# undef _igC_
#elif defined (_lgC_)
# define gCMatrix_colSums lgCMatrix_colSums_i
# define _INT_ans
# define _has_x_slot_
/*Future? # define _has_x_l_slot_ */
# undef _lgC_
#elif defined (_lgC_mn)
# define gCMatrix_colSums lgCMatrix_colSums_d
# define _DOUBLE_ans
# define _has_x_slot_
/*Future? # define _has_x_l_slot_ */
# undef _lgC_mn
#elif defined (_ngC_)
# define gCMatrix_colSums ngCMatrix_colSums_i
# define _INT_ans
/* withOUT 'x' slot */
# undef _ngC_
#elif defined (_ngC_mn)
# define gCMatrix_colSums ngCMatrix_colSums_d
# define _DOUBLE_ans
/* withOUT 'x' slot */
# undef _ngC_mn
#elif defined (_zgC_)
# error "zgC* not yet implemented"
#else
# error "no valid _[dilnz]gC_ option"
#endif
/* - - - - - - - - - - - - - - - - - - - - */
/* Most of this is maybe for the future,
* when cholmod has integer 'x' slot :*/
#ifdef _has_x_d_slot_
# define Type_x_ double
# define STYP_x_ REAL
# define _has_x_slot_
# undef _has_x_d_slot_
#elif defined (_has_x_i_slot_)
# define Type_x_ int
# define STYP_x_ INTEGER
# define _has_x_slot_
# undef _has_x_i_slot_
#elif defined (_has_x_l_slot_)
# define Type_x_ int
# define STYP_x_ LOGICAL
# define _has_x_slot_
# undef _has_x_l_slot_
#endif
/* - - - - - - - - - - - - - - - - - - - - */
#ifdef _DOUBLE_ans
# define SparseResult_class "dsparseVector"
# define Type_ans double
# define STYP_ans REAL
# define NA_ans NA_REAL
# define SXP_ans REALSXP
# define COERCED(x) (x)
#undef _DOUBLE_ans
#elif defined (_INT_ans)
# define SparseResult_class "isparseVector"
# define Type_ans int
# define STYP_ans INTEGER
# define NA_ans NA_INTEGER
# define SXP_ans INTSXP
# define COERCED(x) (Type_ans)(x != 0)
#undef _INT_ans
#else
# error "invalid macro logic"
#endif
/* - - - - - - - - - - - - - - - - - - - - */
#ifdef _has_x_slot_
/* currently have x slot always double (cholmod restriction): */
# define is_NA_x_(u) ISNAN(u)
# define ColSUM_column(_i1_,_i2_,_SUM_) \
if(mn) dnm = cx->nrow; /* denominator for means */ \
for(i = _i1_, _SUM_ = 0; i < _i2_; i++) { \
if (is_NA_x_(xx[i])) { \
if(!na_rm) { \
_SUM_ = NA_ans; \
break; \
} \
/* else: na_rm : skip NAs , */ \
if(mn) /* but decrement denominator */ \
dnm--; \
} else _SUM_ += COERCED(xx[i]); \
} \
if(mn) _SUM_ = (dnm > 0) ? _SUM_/dnm : NA_ans
#else /* no 'x' slot -> no NAs ... */
# define ColSUM_column(_i1_,_i2_,_SUM_) \
_SUM_ = _i2_ - _i1_; \
if(mn) _SUM_ /= cx->nrow
#endif
/* Now the template which depends on the above macros : */
SEXP gCMatrix_colSums(SEXP x, SEXP NArm, SEXP spRes, SEXP trans, SEXP means)
{
int mn = asLogical(means), sp = asLogical(spRes), tr = asLogical(trans);
/* cholmod_sparse: drawback of coercing lgC to double: */
CHM_SP cx = AS_CHM_SP(x);
R_CheckStack();
if (tr) {
cholmod_sparse *cxt = cholmod_transpose(cx, (int)cx->xtype, &c);
cx = cxt;
}
/* everything else *after* the above potential transpose : */
/* Don't declarations here require the C99 standard? Can we assume C99? */
int j, nc = cx->ncol;
int *xp = (int *)(cx -> p);
#ifdef _has_x_slot_
int na_rm = asLogical(NArm), i, dnm = 0/*Wall*/;
double *xx = (double *)(cx -> x);
#endif
SEXP ans = PROTECT(sp ? NEW_OBJECT(MAKE_CLASS(SparseResult_class))
: allocVector(SXP_ans, nc));
if (sp) { /* sparseResult - never allocating length-nc ... */
int nza, i1, i2, p, *ai;
Type_ans *ax;
for (j = 0, nza = 0; j < nc; j++)
if(xp[j] < xp[j + 1])
nza++;
ai = INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nza));
ax = STYP_ans(ALLOC_SLOT(ans, Matrix_xSym, SXP_ans, nza));
SET_SLOT(ans, Matrix_lengthSym, ScalarInteger(nc));
i2 = xp[0];
for (j = 1, p = 0; j <= nc; j++) {
/* j' =j+1, since 'i' slot will be 1-based */
i1 = i2; i2 = xp[j];
if(i1 < i2) {
Type_ans sum;
ColSUM_column(i1,i2, sum);
ai[p] = j;
ax[p++] = sum;
}
}
}
else { /* "numeric" (non sparse) result */
Type_ans *a = STYP_ans(ans);
for (j = 0; j < nc; j++) {
ColSUM_column(xp[j], xp[j + 1], a[j]);
}
}
if (tr) cholmod_free_sparse(&cx, &c);
UNPROTECT(1);
return ans;
}
#undef ColSUM_column
#undef NA_ans
#undef STYP_ans
#undef SXP_ans
#undef SparseResult_class
#undef Type_ans
#undef COERCED
#ifdef _has_x_slot_
# undef NA_x_
# undef Type_x_
# undef STYP_x_
# undef _has_x_slot_
#endif
#undef gCMatrix_colSums