/*------ 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