Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
t_Csparse_subassign.c
Go to the documentation of this file.
1 /*------ Definition of a template for [diln]Csparse_subassign(...) : *
2  * -------- ~~~~~~~~~~~~~~~~~~~~~~
3  * i.e., included several times from ./Csparse.c
4  * ~~~~~~~~~~~
5  *
6  _slot_kind : use the integer codes matching x_slot_kind in ./Mutils.h
7  * ~~~~~~~~
8  */
9 
10 #ifdef _d_Csp_
11 
12 # define Csparse_subassign dCsparse_subassign
13 # define x_CLASSES "dgCMatrix",/* 0 */ "dtCMatrix" /* 1 */
14 # define sparseVECTOR "dsparseVector"
15 # define slot_kind_x 0
16 # define _DOUBLE_x
17 # define _has_x_slot_
18 
19 # undef _d_Csp_
20 
21 #elif defined (_l_Csp_)
22 
23 # define Csparse_subassign lCsparse_subassign
24 # define x_CLASSES "lgCMatrix",/* 0 */ "ltCMatrix" /* 1 */
25 # define sparseVECTOR "lsparseVector"
26 # define slot_kind_x 1
27 # define _LGL_x
28 # define _has_x_slot_
29 
30 # undef _l_Csp_
31 
32 #elif defined (_i_Csp_)
33 
34 # define Csparse_subassign iCsparse_subassign
35 # define x_CLASSES "igCMatrix",/* 0 */ "itCMatrix" /* 1 */
36 # define sparseVECTOR "isparseVector"
37 # define slot_kind_x 2
38 # define _INT_x
39 # define _has_x_slot_
40 
41 # undef _i_Csp_
42 
43 #elif defined (_n_Csp_)
44 
45 # define Csparse_subassign nCsparse_subassign
46 # define x_CLASSES "ngCMatrix",/* 0 */ "ntCMatrix" /* 1 */
47 # define sparseVECTOR "nsparseVector"
48 # define slot_kind_x -1
49 # define _INT_x
50 /* withOUT 'x' slot -- CARE! we assume that this is the *ONLY* case w/o x slot */
51 
52 # undef _n_Csp_
53 
54 #elif defined (_z_Csp_)
55 
56 // # error "zgC* not yet implemented"
57 # define Csparse_subassign zCsparse_subassign
58 # define x_CLASSES "zgCMatrix",/* 0 */ "ztCMatrix" /* 1 */
59 # define sparseVECTOR "zsparseVector"
60 # define slot_kind_x 3
61 # define _CPLX_x
62 # define _has_x_slot_
63 
64 # undef _z_Csp_
65 
66 #else
67 
68 # error "no valid _[dilnz]gC_ option"
69 
70 #endif
71 // -------------------------------------------------
72 
73 #ifdef _DOUBLE_x
74 
75 # define Type_x double
76 # define Type_x_0_init(_VAR_) double _VAR_ = 0.
77 # define Type_x_1_init(_VAR_) double _VAR_ = 1.
78 # define STYP_x REAL
79 # define SXP_x REALSXP
80 #undef _DOUBLE_x
81 
82 #elif defined (_LGL_x)
83 
84 # define Type_x int
85 # define Type_x_0_init(_VAR_) int _VAR_ = 0
86 # define Type_x_1_init(_VAR_) int _VAR_ = 1
87 # define STYP_x LOGICAL
88 # define SXP_x LGLSXP
89 #undef _LGL_x
90 
91 #elif defined (_INT_x)
92 
93 # define Type_x int
94 # define Type_x_0_init(_VAR_) int _VAR_ = 0
95 # define Type_x_1_init(_VAR_) int _VAR_ = 1
96 # define STYP_x INTEGER
97 # define SXP_x INTSXP
98 #undef _INT_x
99 
100 #elif defined (_CPLX_x)
101 
102 # define Type_x Rcomplex
103 # define Type_x_0_init(_VAR_) Rcomplex _VAR_; _VAR_.r = _VAR_.i = 0.
104 # define Type_x_1_init(_VAR_) Rcomplex _VAR_; _VAR_.r = 1.; _VAR_.i = 0.
105 # define STYP_x COMPLEX
106 # define SXP_x CPLXSXP
107 
108 #else
109 # error "invalid macro logic"
110 #endif
111 
112 
113 
124 SEXP Csparse_subassign(SEXP x, SEXP i_, SEXP j_, SEXP value)
125 {
126  // TODO: for other classes consider using a trick as RallocedReal() in ./chm_common.c
127  static const char
128  *valid_cM [] = { // the only ones, for "the moment". FIXME: extend (!)
129  x_CLASSES,
130  ""},
131  // value: assume a "dsparseVector" for now -- slots: (i, length, x)
132  *valid_spv[] = { sparseVECTOR, // = "the one with the same slot-class"
133  // all others: ctype_v slot_kind
134  "nsparseVector",// 1 -1
135  "lsparseVector",// 2 1
136  "isparseVector",// 3 2
137  "dsparseVector",// 4 0
138  "zsparseVector",// 5 3
139  ""};
140 
141  int ctype_x = R_check_class_etc(x, valid_cM),
142  ctype_v = R_check_class_etc(value, valid_spv);
143  if (ctype_x < 0)
144  error(_("invalid class of 'x' in Csparse_subassign()"));
145  if (ctype_v < 0)
146  error(_("invalid class of 'value' in Csparse_subassign()"));
147  Rboolean value_is_nsp = ctype_v == 1;
148 #ifndef _has_x_slot_ // i.e. "n.CMatrix" : sparseVECTOR == "nsparseVector"
149  if(!value_is_nsp) value_is_nsp = (ctype_v == 0);
150 #endif
151 
152  SEXP
153  islot = GET_SLOT(x, Matrix_iSym),
154  dimslot = GET_SLOT(x, Matrix_DimSym),
155  i_cp = PROTECT(coerceVector(i_, INTSXP)),
156  j_cp = PROTECT(coerceVector(j_, INTSXP));
157  // for d.CMatrix and l.CMatrix but not n.CMatrix:
158 
159  int *dims = INTEGER(dimslot),
160  ncol = dims[1], /* nrow = dims[0], */
161  *i = INTEGER(i_cp), len_i = LENGTH(i_cp),
162  *j = INTEGER(j_cp), len_j = LENGTH(j_cp),
163  k,
164  nnz_x = LENGTH(islot);
165  int nnz = nnz_x;
166 
167 #define MATRIX_SUBASSIGN_VERBOSE
168 // Temporary hack for debugging --- remove eventually -- FIXME
169 #ifdef MATRIX_SUBASSIGN_VERBOSE
170  Rboolean verbose = i[0] < 0;
171  if(verbose) {
172  i[0] = -i[0];
173  REprintf("Csparse_subassign() x[i,j] <- val; x is \"%s\"; value \"%s\" is_nsp=%d\n",
174  valid_cM[ctype_x], valid_spv[ctype_v], (int)value_is_nsp);
175  }
176 #endif
177 
178  SEXP val_i_slot, val_x_slot;
179  val_i_slot = PROTECT(coerceVector(GET_SLOT(value, Matrix_iSym), REALSXP));
180  double *val_i = REAL(val_i_slot);
181  int nnz_val = LENGTH(GET_SLOT(value, Matrix_iSym)), n_prot = 4;
182  Type_x *val_x = NULL;
183  if(!value_is_nsp) {
184  if(ctype_v) { // matrix 'x' and 'value' are of different kinds
185  switch((enum x_slot_kind) slot_kind_x) {
186  case x_pattern:// "n"
187  case x_logical:// "l"
188  if(ctype_v >= 3)
189  warning(_("x[] <- val: val is coerced to logical for \"%s\" x"),
190  valid_cM[ctype_x]);
191  break;
192  case x_integer:
193  if(ctype_v >= 4)
194  error(_("x[] <- val: val should be integer or logical, is coerced to integer, for \"%s\" x"),
195  valid_cM[ctype_x]);
196  break;
197  case x_double:
198  case x_complex: // coercion should be tried (and fail for complex -> double) below
199  break;
200  default:
201  error(_("programming error in Csparse_subassign() should never happen"));
202  }
203  // otherwise: "coerce" : as(., <sparseVector>) :
204  val_x_slot = PROTECT(coerceVector(GET_SLOT(value, Matrix_xSym), SXP_x)); n_prot++;
205  val_x = STYP_x(val_x_slot);
206  } else {
207  val_x = STYP_x( GET_SLOT(value, Matrix_xSym));
208  }
209  }
210  int64_t len_val = (int64_t) asReal(GET_SLOT(value, Matrix_lengthSym));
211  /* llen_i = (int64_t) len_i; */
212 
213  SEXP ans;
214  /* Instead of simple "duplicate": PROTECT(ans = duplicate(x)) , build up: */
215  // Assuming that ans will have the same basic Matrix type as x :
216  ans = PROTECT(NEW_OBJECT(MAKE_CLASS(valid_cM[ctype_x])));
217  SET_SLOT(ans, Matrix_DimSym, duplicate(dimslot));
218  slot_dup(ans, x, Matrix_DimNamesSym);
219  slot_dup(ans, x, Matrix_pSym);
220  SEXP r_pslot = GET_SLOT(ans, Matrix_pSym);
221  // and assign the i- and x- slots at the end, as they are potentially modified
222  // not just in content, but also in their *length*
223  int *rp = INTEGER(r_pslot),
224  *ri = Calloc(nnz_x, int); // to contain the final i - slot
225  Memcpy(ri, INTEGER(islot), nnz_x);
226  Type_x_0_init(z_ans);
227  Type_x_1_init(one_ans);
228 #ifdef _has_x_slot_
229  Type_x *rx = Calloc(nnz_x, Type_x); // to contain the final x - slot
230  Memcpy(rx, STYP_x(GET_SLOT(x, Matrix_xSym)), nnz_x);
231 #endif
232  // NB: nnz_x : will always be the "current allocated length" of (i, x) slots
233  // -- nnz : the current *used* length; always nnz <= nnz_x
234 
235  int jj, j_val = 0; // in "running" conceptionally through all value[i+ jj*len_i]
236  // values, we are "below"/"before" the (j_val)-th non-zero one.
237  // e.g. if value = (0,0,...,0), have nnz_val == 0, j_val must remain == 0
238  int64_t ii_val;// == "running" index (i + jj*len_i) % len_val for value[]
239  for(jj = 0, ii_val=0; jj < len_j; jj++) {
240  int j__ = j[jj];
241  /* int64_t j_l = jj * llen_i; */
242  R_CheckUserInterrupt();
243  for(int ii = 0; ii < len_i; ii++, ii_val++) {
244  int i__ = i[ii], p1, p2;
245  if(nnz_val && ii_val >= len_val) { // "recycle" indexing into value[]
246  ii_val -= len_val; // = (ii + jj*len_i) % len_val
247  j_val = 0;
248  }
249  int64_t ii_v1;//= ii_val + 1;
250  Type_x v, /* := value[(ii + j_l) % len_val]
251  = .sparseVector_sub((ii + j_l) % len_val,
252  nnz_val, val_i, val_x, len_val)
253  */
254  M_ij;
255  int ind;
256  Rboolean have_entry = FALSE;
257 
258  // note that rp[]'s may have *changed* even when 'j' remained!
259  // "FIXME": do this only *when* rp[] has changed
260  p1 = rp[j__], p2 = rp[j__ + 1];
261 
262  // v := value[(ii + j_l) % len_val] = value[ii_val]
263  v = z_ans;
264  if(j_val < nnz_val) { // maybe find v := non-zero value[ii_val]
265  ii_v1 = ii_val + 1;
266  if(ii_v1 < val_i[j_val]) { // typical case: are still in zero-stretch
267  // v = z_ans (== 0)
268  } else if(ii_v1 == val_i[j_val]) { // have a match
269  v = (value_is_nsp) ? one_ans : val_x[j_val];
270  j_val++;// from now on, look at the next non-zero entry
271  } else { // ii_v1 > val_i[j_val]
272  REprintf("programming thinko in Csparse_subassign(*, i=%d,j=%d): ii_v=%d, v@i[j_val=%ld]=%g\n",
273  i__,j__, ii_v1, j_val, val_i[j_val]);
274  j_val++;// from now on, look at the next non-zero entry
275  }
276  }
277  // --------------- M_ij := getM(i., j.) --------------------------------
278  M_ij = z_ans; // as in ./t_sparseVector.c
279  for(ind = p1; ind < p2; ind++) {
280  if(ri[ind] >= i__) {
281  if(ri[ind] == i__) {
282 #ifdef _has_x_slot_
283  M_ij = rx[ind];
284 #else
285  M_ij = 1;
286 #endif
287 #ifdef MATRIX_SUBASSIGN_VERBOSE
288  if(verbose)
289  REprintf("have entry x[%d, %d] = %g\n", i__, j__,
290 # ifdef _CPLX_x
291  (double)M_ij.r);
292 # else
293  (double)M_ij);
294 # endif
295 #endif
296  have_entry = TRUE;
297  } else { // ri[ind] > i__
298 #ifdef MATRIX_SUBASSIGN_VERBOSE
299  if(verbose)
300  REprintf("@i > i__ = %d --> ind-- = %d\n", i__, ind);
301 #endif
302  }
303  break;
304  }
305  }
306 
307  //-- R: if(getM(i., j.) != (v <- getV(ii, jj)))
308 
309  // if(contents differ) ==> value needs to be changed :
310 #ifdef _CPLX_x
311  if(M_ij.r != v.r || M_ij.i != v.i) {
312 #else
313  if(M_ij != v) {
314 #endif
315 
316 #ifdef MATRIX_SUBASSIGN_VERBOSE
317  if(verbose)
318  REprintf("setting x[%d, %d] <- %g", i__,j__,
319 # ifdef _CPLX_x
320  (double) v.r);
321 # else
322  (double) v);
323 # endif
324 #endif
325  // (otherwise: nothing to do):
326  // setM(i__, j__, v)
327  // ----------------------------------------------------------
328 
329 #ifndef _has_x_slot_
330  if(v == z_ans) {
331  // Case I ----- remove x[i, j] = M_ij which we know is *non*-zero
332 // BUT it is more efficient (memory-management!) *NOT* to remove,
334 // currently using drop0() in R code
335  // we know : have_entry = TRUE ;
336  // ri[ind] == i__; M_ij = rx[ind];
337 #ifdef MATRIX_SUBASSIGN_VERBOSE
338  if(verbose)
339  REprintf(" rm ind=%d\n", ind);
340 #endif
341  // remove the 'ind'-th element from x@i and x@x :
342  nnz-- ;
343  for(k=ind; k < nnz; k++) {
344  ri[k] = ri[k+1];
345 #ifdef _has_x_slot_
346  rx[k] = rx[k+1];
347 #endif
348  }
349  for(k=j__ + 1; k <= ncol; k++) {
350  rp[k] = rp[k] - 1;
351  }
352  }
353  else
354 #endif
355  if(have_entry) {
356  // Case II ----- replace (non-empty) x[i,j] by v -------
357 #ifdef MATRIX_SUBASSIGN_VERBOSE
358  if(verbose)
359  REprintf(" repl. ind=%d\n", ind);
360 #endif
361 #ifdef _has_x_slot_
362  rx[ind] = v;
363 #endif
364  } else {
365  // Case III ----- v != 0 : insert v into "empty" x[i,j] ----
366 
367  // extend the i and x slot by one entry : ---------------------
368  if(nnz+1 > nnz_x) { // need to reallocate:
369 #ifdef MATRIX_SUBASSIGN_VERBOSE
370  if(verbose) REprintf(" Realloc()ing: nnz_x=%d", nnz_x);
371 #endif
372  // do it "only" 1x,..4x at the very most increasing by the
373  // nnz-length of "value":
374  nnz_x += (1 + nnz_val / 4);
375 #ifdef MATRIX_SUBASSIGN_VERBOSE
376  if(verbose) REprintf("(nnz_v=%d) --> %d ", nnz_val, nnz_x);
377 #endif
378  // C doc on realloc() says that the old content is *preserve*d
379  ri = Realloc(ri, nnz_x, int);
380 #ifdef _has_x_slot_
381  rx = Realloc(rx, nnz_x, Type_x);
382 #endif
383  }
384  // 3) fill them ...
385 
386  int i1 = ind;
387 #ifdef MATRIX_SUBASSIGN_VERBOSE
388  if(verbose)
389  REprintf(" INSERT p12=(%d,%d) -> ind=%d -> i1 = %d\n",
390  p1,p2, ind, i1);
391 #endif
392 
393  // shift the "upper values" *before* the insertion:
394  for(int l = nnz-1; l >= i1; l--) {
395  ri[l+1] = ri[l];
396 #ifdef _has_x_slot_
397  rx[l+1] = rx[l];
398 #endif
399  }
400  ri[i1] = i__;
401 #ifdef _has_x_slot_
402  rx[i1] = v;
403 #endif
404  nnz++;
405  // the columns j "right" of the current one :
406  for(k=j__ + 1; k <= ncol; k++)
407  rp[k]++;
408  }
409  }
410 #ifdef MATRIX_SUBASSIGN_VERBOSE
411  else if(verbose) REprintf("M_ij == v = %g\n",
412 # ifdef _CPLX_x
413  (double) v.r);
414 # else
415  (double) v);
416 # endif
417 #endif
418  }// for( ii )
419  }// for( jj )
420 
421  if(ctype_x == 1) { // triangularMatrix: copy the 'diag' and 'uplo' slots
422  slot_dup(ans, x, Matrix_uploSym);
423  slot_dup(ans, x, Matrix_diagSym);
424  }
425  // now assign the i- and x- slots, free memory and return :
426  Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)), ri, nnz);
427 #ifdef _has_x_slot_
428  Memcpy( STYP_x(ALLOC_SLOT(ans, Matrix_xSym, SXP_x, nnz)), rx, nnz);
429  Free(rx);
430 #endif
431  Free(ri);
432  UNPROTECT(n_prot);
433  return ans;
434 }
435 
436 #undef Csparse_subassign
437 #undef x_CLASSES
438 #undef sparseVECTOR
439 
440 #undef Type_x
441 #undef STYP_x
442 #undef SXP_x
443 #undef Type_x_0_init
444 #undef Type_x_1_init
445 #undef _has_x_slot_
446 #undef slot_kind_x
447 
448 #ifdef _CPLX_x
449 # undef _CPLX_x
450 #endif
451 
452 
453 
SEXP Matrix_DimSym
Definition: Syms.h:2
x_slot_kind
Definition: Mutils.h:184
SEXP Matrix_xSym
Definition: Syms.h:2
#define slot_dup(dest, src, sym)
Definition: Mutils.h:149
SEXP Matrix_lengthSym
Definition: Syms.h:2
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
SEXP Matrix_uploSym
Definition: Syms.h:2
#define _(String)
Definition: Mutils.h:32
SEXP Matrix_iSym
Definition: Syms.h:2
SEXP Csparse_subassign(SEXP x, SEXP i_, SEXP j_, SEXP value)
Subassignment: x[i,j] <- value.
SEXP Matrix_pSym
Definition: Syms.h:2
SEXP Matrix_diagSym
Definition: Syms.h:2
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