Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
t_gCMatrix_colSums.c
Go to the documentation of this file.
1 /*------ Definition of a template for [diln]gCMatrix_colsums(...) : *
2  * -------- ~~~~~~~~~~~~~~~~~~~~~~
3  * i.e., included several times from ./dgCMatrix.c
4  * ~~~~~~~~~~~~~
5  */
6 
7 
8 /* for all cases with an 'x' slot -- i.e. almost all cases ;
9  * just redefine this in the other cases:
10  */
11 
12 #ifdef _dgC_
13 
14 # define gCMatrix_colSums dgCMatrix_colSums
15 # define _DOUBLE_ans
16 # define _has_x_slot_
17 /*Future? # define _has_x_d_slot_ */
18 # undef _dgC_
19 
20 #elif defined (_igC_)
21 
22 # define gCMatrix_colSums igCMatrix_colSums
23 # define _DOUBLE_ans
24 # define _has_x_slot_
25 /*Future? # define _has_x_d_slot_ */
26 # undef _igC_
27 
28 #elif defined (_lgC_)
29 
30 # define gCMatrix_colSums lgCMatrix_colSums_i
31 # define _INT_ans
32 # define _has_x_slot_
33 /*Future? # define _has_x_l_slot_ */
34 # undef _lgC_
35 
36 #elif defined (_lgC_mn)
37 
38 # define gCMatrix_colSums lgCMatrix_colSums_d
39 # define _DOUBLE_ans
40 # define _has_x_slot_
41 /*Future? # define _has_x_l_slot_ */
42 # undef _lgC_mn
43 
44 #elif defined (_ngC_)
45 
46 # define gCMatrix_colSums ngCMatrix_colSums_i
47 # define _INT_ans
48  /* withOUT 'x' slot */
49 # undef _ngC_
50 
51 #elif defined (_ngC_mn)
52 
53 # define gCMatrix_colSums ngCMatrix_colSums_d
54 # define _DOUBLE_ans
55  /* withOUT 'x' slot */
56 # undef _ngC_mn
57 
58 #elif defined (_zgC_)
59 
60 # error "zgC* not yet implemented"
61 
62 #else
63 
64 # error "no valid _[dilnz]gC_ option"
65 
66 #endif
67 
68 /* - - - - - - - - - - - - - - - - - - - - */
69 
70 /* Most of this is maybe for the future,
71  * when cholmod has integer 'x' slot :*/
72 #ifdef _has_x_d_slot_
73 
74 # define Type_x_ double
75 # define STYP_x_ REAL
76 # define _has_x_slot_
77 # undef _has_x_d_slot_
78 
79 #elif defined (_has_x_i_slot_)
80 
81 # define Type_x_ int
82 # define STYP_x_ INTEGER
83 # define _has_x_slot_
84 # undef _has_x_i_slot_
85 
86 #elif defined (_has_x_l_slot_)
87 
88 # define Type_x_ int
89 # define STYP_x_ LOGICAL
90 # define _has_x_slot_
91 # undef _has_x_l_slot_
92 
93 #endif
94 
95 /* - - - - - - - - - - - - - - - - - - - - */
96 
97 #ifdef _DOUBLE_ans
98 
99 # define SparseResult_class "dsparseVector"
100 # define Type_ans double
101 # define STYP_ans REAL
102 # define NA_ans NA_REAL
103 # define SXP_ans REALSXP
104 # define COERCED(x) (x)
105 #undef _DOUBLE_ans
106 
107 #elif defined (_INT_ans)
108 
109 # define SparseResult_class "isparseVector"
110 # define Type_ans int
111 # define STYP_ans INTEGER
112 # define NA_ans NA_INTEGER
113 # define SXP_ans INTSXP
114 # define COERCED(x) (Type_ans)(x != 0)
115 #undef _INT_ans
116 
117 #else
118 # error "invalid macro logic"
119 #endif
120 
121 /* - - - - - - - - - - - - - - - - - - - - */
122 
123 #ifdef _has_x_slot_
124 
125 /* currently have x slot always double (cholmod restriction): */
126 # define is_NA_x_(u) ISNAN(u)
127 
128 # define ColSUM_column(_i1_,_i2_,_SUM_) \
129  if(mn) dnm = cx->nrow; /* denominator for means */ \
130  for(i = _i1_, _SUM_ = 0; i < _i2_; i++) { \
131  if (is_NA_x_(xx[i])) { \
132  if(!na_rm) { \
133  _SUM_ = NA_ans; \
134  break; \
135  } \
136  /* else: na_rm : skip NAs , */ \
137  if(mn) /* but decrement denominator */ \
138  dnm--; \
139  } else _SUM_ += COERCED(xx[i]); \
140  } \
141  if(mn) _SUM_ = (dnm > 0) ? _SUM_/dnm : NA_ans
142 
143 #else /* no 'x' slot -> no NAs ... */
144 
145 # define ColSUM_column(_i1_,_i2_,_SUM_) \
146  _SUM_ = _i2_ - _i1_; \
147  if(mn) _SUM_ /= cx->nrow
148 #endif
149 
150 /* Now the template which depends on the above macros : */
151 
160 SEXP gCMatrix_colSums(SEXP x, SEXP NArm, SEXP spRes, SEXP trans, SEXP means)
161 {
162  int mn = asLogical(means), sp = asLogical(spRes), tr = asLogical(trans);
163  /* cholmod_sparse: drawback of coercing lgC to double: */
164  CHM_SP cx = AS_CHM_SP__(x);
165  R_CheckStack();
166 
167  if (tr) {
168  cholmod_sparse *cxt = cholmod_transpose(cx, (int)cx->xtype, &c);
169  cx = cxt;
170  }
171 
172  /* everything else *after* the above potential transpose : */
173 
174  int j, nc = cx->ncol;
175  int *xp = (int *)(cx -> p);
176 #ifdef _has_x_slot_
177  int na_rm = asLogical(NArm), // can have NAs only with an 'x' slot
178  i, dnm = 0/*Wall*/;
179  double *xx = (double *)(cx -> x);
180 #endif
181  // result value: sparseResult (==> "*sparseVector") or dense (atomic)vector
182  SEXP ans = PROTECT(sp ? NEW_OBJECT(MAKE_CLASS(SparseResult_class))
183  : allocVector(SXP_ans, nc));
184  if (sp) { // sparseResult, i.e. *sparseVector (never allocating length-nc)
185  int nza, i1, i2, p, *ai;
186  Type_ans *ax;
187 
188  for (j = 0, nza = 0; j < nc; j++)
189  if(xp[j] < xp[j + 1])
190  nza++;
191 
192  ai = INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nza));
193  ax = STYP_ans(ALLOC_SLOT(ans, Matrix_xSym, SXP_ans, nza));
194 
195  SET_SLOT(ans, Matrix_lengthSym, ScalarInteger(nc));
196 
197  i2 = xp[0];
198  for (j = 1, p = 0; j <= nc; j++) {
199  /* j' =j+1, since 'i' slot will be 1-based */
200  i1 = i2; i2 = xp[j];
201  if(i1 < i2) {
202  Type_ans sum;
203  ColSUM_column(i1,i2, sum);
204 
205  ai[p] = j;
206  ax[p++] = sum;
207  }
208  }
209  }
210  else { /* "numeric" (non sparse) result */
211  Type_ans *a = STYP_ans(ans);
212  for (j = 0; j < nc; j++) {
213  ColSUM_column(xp[j], xp[j + 1], a[j]);
214  }
215  }
216 
217  if (tr) cholmod_free_sparse(&cx, &c);
218  if (!sp) {
219  SEXP nms = VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), tr ? 0 : 1);
220  if (!isNull(nms))
221  setAttrib(ans, R_NamesSymbol, duplicate(nms));
222  }
223  UNPROTECT(1);
224  return ans;
225 }
226 
227 #undef ColSUM_column
228 
229 #undef NA_ans
230 #undef STYP_ans
231 #undef SXP_ans
232 #undef SparseResult_class
233 #undef Type_ans
234 
235 #undef COERCED
236 
237 #ifdef _has_x_slot_
238 # undef NA_x_
239 # undef Type_x_
240 # undef STYP_x_
241 # undef _has_x_slot_
242 #endif
243 
244 #undef gCMatrix_colSums
#define AS_CHM_SP__(x)
Definition: chm_common.h:49
cholmod_sparse * CHM_SP
Definition: chm_common.h:25
SEXP Matrix_xSym
Definition: Syms.h:2
SEXP Matrix_lengthSym
Definition: Syms.h:2
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
#define ColSUM_column(_i1_, _i2_, _SUM_)
SEXP Matrix_iSym
Definition: Syms.h:2
SEXP gCMatrix_colSums(SEXP x, SEXP NArm, SEXP spRes, SEXP trans, SEXP means)
colSums(), colMeans(), rowSums() and rowMeans() for all sparce *gCMatrix()es
cholmod_common c
Definition: chm_common.c:15
static R_INLINE SEXP ALLOC_SLOT(SEXP obj, SEXP nm, SEXPTYPE type, int length)
Allocate an SEXP of given type and length, assign it as slot nm in the object, and return the SEXP...
Definition: Mutils.h:240