Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
Csparse.c
Go to the documentation of this file.
1 
6 #include "Csparse.h"
7 #include "Tsparse.h"
8 #include "chm_common.h"
9 
11 Rboolean isValid_Csparse(SEXP x)
12 {
13  /* NB: we do *NOT* check a potential 'x' slot here, at all */
14  SEXP pslot = GET_SLOT(x, Matrix_pSym),
15  islot = GET_SLOT(x, Matrix_iSym);
16  int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), j,
17  nrow = dims[0],
18  ncol = dims[1],
19  *xp = INTEGER(pslot),
20  *xi = INTEGER(islot);
21 
22  if (length(pslot) != dims[1] + 1)
23  return FALSE;
24  if (xp[0] != 0)
25  return FALSE;
26  if (length(islot) < xp[ncol]) /* allow larger slots from over-allocation!*/
27  return FALSE;
28  for (j = 0; j < xp[ncol]; j++) {
29  if (xi[j] < 0 || xi[j] >= nrow)
30  return FALSE;
31  }
32  for (j = 0; j < ncol; j++) {
33  if (xp[j] > xp[j + 1])
34  return FALSE;
35  }
36  return TRUE;
37 }
38 
39 SEXP Csparse_validate(SEXP x) {
40  return Csparse_validate_(x, FALSE);
41 }
42 
43 
44 #define _t_Csparse_validate
45 #include "t_Csparse_validate.c"
46 
47 #define _t_Csparse_sort
48 #include "t_Csparse_validate.c"
49 
50 // R: .validateCsparse(x, sort.if.needed = FALSE) :
51 SEXP Csparse_validate2(SEXP x, SEXP maybe_modify) {
52  return Csparse_validate_(x, asLogical(maybe_modify));
53 }
54 
55 // R: Matrix:::.sortCsparse(x) :
56 SEXP Csparse_sort (SEXP x) {
57  int ok = Csparse_sort_2(x, TRUE); // modifying x directly
58  if(!ok) warning(_("Csparse_sort(x): x is not a valid (apart from sorting) CsparseMatrix"));
59  return x;
60 }
61 
62 SEXP Rsparse_validate(SEXP x)
63 {
64  /* NB: we do *NOT* check a potential 'x' slot here, at all */
65  SEXP pslot = GET_SLOT(x, Matrix_pSym),
66  jslot = GET_SLOT(x, Matrix_jSym);
67  Rboolean sorted, strictly;
68  int i, k,
69  *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
70  nrow = dims[0],
71  ncol = dims[1],
72  *xp = INTEGER(pslot),
73  *xj = INTEGER(jslot);
74 
75  if (length(pslot) != dims[0] + 1)
76  return mkString(_("slot p must have length = nrow(.) + 1"));
77  if (xp[0] != 0)
78  return mkString(_("first element of slot p must be zero"));
79  if (length(jslot) < xp[nrow]) /* allow larger slots from over-allocation!*/
80  return
81  mkString(_("last element of slot p must match length of slots j and x"));
82  for (i = 0; i < length(jslot); i++) {
83  if (xj[i] < 0 || xj[i] >= ncol)
84  return mkString(_("all column indices must be between 0 and ncol-1"));
85  }
86  sorted = TRUE; strictly = TRUE;
87  for (i = 0; i < nrow; i++) {
88  if (xp[i] > xp[i+1])
89  return mkString(_("slot p must be non-decreasing"));
90  if(sorted)
91  for (k = xp[i] + 1; k < xp[i + 1]; k++) {
92  if (xj[k] < xj[k - 1])
93  sorted = FALSE;
94  else if (xj[k] == xj[k - 1])
95  strictly = FALSE;
96  }
97  }
98  if (!sorted)
99  /* cannot easily use cholmod_sort(.) ... -> "error out" :*/
100  return mkString(_("slot j is not increasing inside a column"));
101  else if(!strictly) /* sorted, but not strictly */
102  return mkString(_("slot j is not *strictly* increasing inside a column"));
103 
104  return ScalarLogical(1);
105 }
106 
121 SEXP Csparse_to_dense(SEXP x, SEXP symm_or_tri)
122 {
123  Rboolean is_sym, is_tri;
124  int is_sym_or_tri = asInteger(symm_or_tri),
125  ctype = 0; // <- default = "dgC"
126  static const char *valid[] = { MATRIX_VALID_Csparse, ""};
127  if(is_sym_or_tri == NA_INTEGER) { // find if is(x, "symmetricMatrix") :
128  ctype = R_check_class_etc(x, valid);
129  is_sym = (ctype % 3 == 1);
130  is_tri = (ctype % 3 == 2);
131  } else {
132  is_sym = is_sym_or_tri > 0;
133  is_tri = is_sym_or_tri < 0;
134  // => both are FALSE iff is_.. == 0
135  if(is_sym || is_tri)
136  ctype = R_check_class_etc(x, valid);
137  }
138  CHM_SP chxs = AS_CHM_SP__(x);// -> chxs->stype = +- 1 <==> symmetric
139  R_CheckStack();
140  if(is_tri && *diag_P(x) == 'U') { // ==> x := diagU2N(x), directly for chxs
141  CHM_SP eye = cholmod_speye(chxs->nrow, chxs->ncol, chxs->xtype, &c);
142  double one[] = {1, 0};
143  CHM_SP ans = cholmod_add(chxs, eye, one, one,
144  /* values: */ ((ctype / 3) != 2), // TRUE iff not "nMatrix"
145  TRUE, &c);
146  cholmod_free_sparse(&eye, &c);
147  chxs = cholmod_copy_sparse(ans, &c);
148  cholmod_free_sparse(&ans, &c);
149  }
150  /* The following loses the symmetry property, since cholmod_dense has none,
151  * BUT, much worse (FIXME!), it also transforms CHOLMOD_PATTERN ("n") matrices
152  * to numeric (CHOLMOD_REAL) ones {and we "revert" via chm_dense_to_SEXP()}: */
153  CHM_DN chxd = cholmod_sparse_to_dense(chxs, &c);
154  /* FIXME: The above FAILS for prod(dim(.)) > INT_MAX
155  /* TODO: use cholmod_l_* but also the 'cl' global ==> many changes in chm_common.[ch]
156  *
157  * >>>>>>>>>>> TODO <<<<<<<<<<<<
158  * CHM_DN chxd = cholmod_l_sparse_to_dense(chxs, &cl); */
159  // ^^^ important when prod(dim(.)) > INT_MAX
160  int Rkind = (chxs->xtype == CHOLMOD_PATTERN)? -1 : Real_kind(x);
161 
162  SEXP ans = chm_dense_to_SEXP(chxd, 1, Rkind, GET_SLOT(x, Matrix_DimNamesSym),
163  /* transp: */ FALSE);
164  // -> a [dln]geMatrix
165  if(is_sym) { // ==> want [dln]syMatrix
166  const char cl1 = class_P(ans)[0];
167  PROTECT(ans);
168  SEXP aa = PROTECT(NEW_OBJECT(MAKE_CLASS((cl1 == 'd') ? "dsyMatrix" :
169  ((cl1 == 'l') ? "lsyMatrix" : "nsyMatrix"))));
170  // No need to duplicate() as slots of ans are freshly allocated and ans will not be used
171  SET_SLOT(aa, Matrix_xSym, GET_SLOT(ans, Matrix_xSym));
172  SET_SLOT(aa, Matrix_DimSym, GET_SLOT(ans, Matrix_DimSym));
173  SET_SLOT(aa, Matrix_DimNamesSym,GET_SLOT(ans, Matrix_DimNamesSym));
174  SET_SLOT(aa, Matrix_uploSym, mkString((chxs->stype > 0) ? "U" : "L"));
175  UNPROTECT(2);
176  return aa;
177  }
178  else if(is_tri) { // ==> want [dln]trMatrix
179  const char cl1 = class_P(ans)[0];
180  PROTECT(ans);
181  SEXP aa = PROTECT(NEW_OBJECT(MAKE_CLASS((cl1 == 'd') ? "dtrMatrix" :
182  ((cl1 == 'l') ? "ltrMatrix" : "ntrMatrix"))));
183  // No need to duplicate() as slots of ans are freshly allocated and ans will not be used
184  SET_SLOT(aa, Matrix_xSym, GET_SLOT(ans, Matrix_xSym));
185  SET_SLOT(aa, Matrix_DimSym, GET_SLOT(ans, Matrix_DimSym));
186  SET_SLOT(aa, Matrix_DimNamesSym,GET_SLOT(ans, Matrix_DimNamesSym));
187  slot_dup(aa, x, Matrix_uploSym);
188  /* already by NEW_OBJECT(..) above:
189  SET_SLOT(aa, Matrix_diagSym, mkString("N")); */
190  UNPROTECT(2);
191  return aa;
192  }
193  else
194  return ans;
195 }
196 
197 // FIXME: do not go via CHM (should not be too hard, to just *drop* the x-slot, right?
198 SEXP Csparse2nz(SEXP x, Rboolean tri)
199 {
200  CHM_SP chxs = AS_CHM_SP__(x);
201  CHM_SP chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c);
202  R_CheckStack();
203 
204  return chm_sparse_to_SEXP(chxcp, 1/*do_free*/,
205  tri ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
206  /* Rkind: pattern */ 0,
207  /* diag = */ tri ? diag_P(x) : "",
208  GET_SLOT(x, Matrix_DimNamesSym));
209 }
210 SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri)
211 {
212  int tr_ = asLogical(tri);
213  if(tr_ == NA_LOGICAL) {
214  warning(_("Csparse_to_nz_pattern(x, tri = NA): 'tri' is taken as TRUE"));
215  tr_ = TRUE;
216  }
217  return Csparse2nz(x, (Rboolean) tr_);
218 }
219 
220 // n.CMatrix --> [dli].CMatrix (not going through CHM!)
221 SEXP nz_pattern_to_Csparse(SEXP x, SEXP res_kind)
222 {
223  return nz2Csparse(x, asInteger(res_kind));
224 }
225 
226 // n.CMatrix --> [dli].CMatrix (not going through CHM!)
227 // NOTE: use chm_MOD_xtype(() to change type of 'cholmod_sparse' matrix
228 SEXP nz2Csparse(SEXP x, enum x_slot_kind r_kind)
229 {
230  const char *cl_x = class_P(x);
231  // quick check - if ok, fast
232  if(cl_x[0] != 'n' || cl_x[2] != 'C') {
233  // e.g. class = "A", from setClass("A", contains = "ngCMatrix")
234  static const char *valid[] = { MATRIX_VALID_nCsparse, ""};
235  int ctype = R_check_class_etc(x, valid);
236  if(ctype < 0)
237  error(_("not a 'n.CMatrix'"));
238  else // fine : get a valid cl_x class_P()-like string :
239  cl_x = valid[ctype];
240  }
241  int nnz = LENGTH(GET_SLOT(x, Matrix_iSym));
242  SEXP ans;
243  char *ncl = alloca(strlen(cl_x) + 1); /* not much memory required */
244  strcpy(ncl, cl_x);
245  double *dx_x; int *ix_x;
246  ncl[0] = (r_kind == x_double ? 'd' :
247  (r_kind == x_logical ? 'l' :
248  /* else (for now): r_kind == x_integer : */ 'i'));
249  PROTECT(ans = NEW_OBJECT(MAKE_CLASS(ncl)));
250  // create a correct 'x' slot:
251  switch(r_kind) {
252  int i;
253  case x_double: // 'd'
254  dx_x = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz));
255  for (i=0; i < nnz; i++) dx_x[i] = 1.;
256  break;
257  case x_logical: // 'l'
258  ix_x = LOGICAL(ALLOC_SLOT(ans, Matrix_xSym, LGLSXP, nnz));
259  for (i=0; i < nnz; i++) ix_x[i] = TRUE;
260  break;
261  case x_integer: // 'i'
262  ix_x = INTEGER(ALLOC_SLOT(ans, Matrix_xSym, INTSXP, nnz));
263  for (i=0; i < nnz; i++) ix_x[i] = 1;
264  break;
265 
266  default:
267  error(_("nz2Csparse(): invalid/non-implemented r_kind = %d"),
268  r_kind);
269  }
270 
271  // now copy all other slots :
272  slot_dup(ans, x, Matrix_iSym);
273  slot_dup(ans, x, Matrix_pSym);
274  slot_dup(ans, x, Matrix_DimSym);
275  slot_dup(ans, x, Matrix_DimNamesSym);
276  if(ncl[1] != 'g') { // symmetric or triangular ...
279  }
280  UNPROTECT(1);
281  return ans;
282 }
283 
284 SEXP Csparse_to_matrix(SEXP x, SEXP chk, SEXP symm)
285 {
286  int is_sym = asLogical(symm);
287  if(is_sym == NA_LOGICAL) { // find if is(x, "symmetricMatrix") :
288  static const char *valid[] = { MATRIX_VALID_Csparse, ""};
289  int ctype = R_check_class_etc(x, valid);
290  is_sym = (ctype % 3 == 1);
291  }
292  return chm_dense_to_matrix(
293  cholmod_sparse_to_dense(AS_CHM_SP2(x, asLogical(chk)), &c),
294  1 /*do_free*/,
295  (is_sym
296  ? symmetric_DimNames(GET_SLOT(x, Matrix_DimNamesSym))
297  : GET_SLOT(x, Matrix_DimNamesSym)));
298 }
299 
300 SEXP Csparse_to_vector(SEXP x)
301 {
302  return chm_dense_to_vector(cholmod_sparse_to_dense(AS_CHM_SP__(x), &c), 1);
303 }
304 
305 SEXP Csparse_to_Tsparse(SEXP x, SEXP tri)
306 {
307  CHM_SP chxs = AS_CHM_SP__(x);
308  CHM_TR chxt = cholmod_sparse_to_triplet(chxs, &c);
309  int tr = asLogical(tri);
310  int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
311  R_CheckStack();
312 
313  return chm_triplet_to_SEXP(chxt, 1,
314  tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
315  Rkind, tr ? diag_P(x) : "",
316  GET_SLOT(x, Matrix_DimNamesSym));
317 }
318 
319 SEXP Csparse_to_tCsparse(SEXP x, SEXP uplo, SEXP diag)
320 {
321  CHM_SP chxs = AS_CHM_SP__(x);
322  int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
323  R_CheckStack();
324  return chm_sparse_to_SEXP(chxs, /* dofree = */ 0,
325  /* uploT = */ (*CHAR(asChar(uplo)) == 'U')? 1: -1,
326  Rkind, /* diag = */ CHAR(STRING_ELT(diag, 0)),
327  GET_SLOT(x, Matrix_DimNamesSym));
328 }
329 
330 SEXP Csparse_to_tTsparse(SEXP x, SEXP uplo, SEXP diag)
331 {
332  CHM_SP chxs = AS_CHM_SP__(x);
333  CHM_TR chxt = cholmod_sparse_to_triplet(chxs, &c);
334  int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
335  R_CheckStack();
336  return chm_triplet_to_SEXP(chxt, 1,
337  /* uploT = */ (*CHAR(asChar(uplo)) == 'U')? 1: -1,
338  Rkind, /* diag = */ CHAR(STRING_ELT(diag, 0)),
339  GET_SLOT(x, Matrix_DimNamesSym));
340 }
341 
342 
344 {
345  CHM_SP chx = AS_CHM_SP__(x), chgx;
346  int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
347  R_CheckStack();
348 
349  if (!(chx->stype))
350  error(_("Nonsymmetric matrix in Csparse_symmetric_to_general"));
351  chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);
352  /* xtype: pattern, "real", complex or .. */
353  return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",
354  symmetric_DimNames(GET_SLOT(x, Matrix_DimNamesSym)));
355 }
356 
357 SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo, SEXP sym_dmns)
358 {
359  int *adims = INTEGER(GET_SLOT(x, Matrix_DimSym)), n = adims[0];
360  if(n != adims[1]) {
361  error(_("Csparse_general_to_symmetric(): matrix is not square!"));
362  return R_NilValue; /* -Wall */
363  }
364  CHM_SP chx = AS_CHM_SP__(x), chgx;
365  int uploT = (*CHAR(asChar(uplo)) == 'U') ? 1 : -1;
366  int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
367  R_CheckStack();
368  chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c);
369 
370  SEXP dns = GET_SLOT(x, Matrix_DimNamesSym);
371  if(asLogical(sym_dmns))
372  dns = symmetric_DimNames(dns);
373  else if((!isNull(VECTOR_ELT(dns, 0)) &&
374  !isNull(VECTOR_ELT(dns, 1))) ||
375  !isNull(getAttrib(dns, R_NamesSymbol))) {
376  /* symmetrize them if both are not NULL
377  * or names(dimnames(.)) is asymmetric : */
378  dns = PROTECT(duplicate(dns));
379  if(!equal_string_vectors(VECTOR_ELT(dns, 0),
380  VECTOR_ELT(dns, 1))) {
381  if(uploT == 1)
382  SET_VECTOR_ELT(dns, 0, VECTOR_ELT(dns,1));
383  else
384  SET_VECTOR_ELT(dns, 1, VECTOR_ELT(dns,0));
385  }
386  SEXP nms_dns = getAttrib(dns, R_NamesSymbol);
387  if(!isNull(nms_dns) && // names(dimnames(.)) :
388  !R_compute_identical(STRING_ELT(nms_dns, 0),
389  STRING_ELT(nms_dns, 1), 16)) {
390  if(uploT == 1)
391  SET_STRING_ELT(nms_dns, 0, STRING_ELT(nms_dns,1));
392  else
393  SET_STRING_ELT(nms_dns, 1, STRING_ELT(nms_dns,0));
394  setAttrib(dns, R_NamesSymbol, nms_dns);
395  }
396  UNPROTECT(1);
397  }
398  /* xtype: pattern, "real", complex or .. */
399  return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", dns);
400 }
401 
402 SEXP Csparse_transpose(SEXP x, SEXP tri)
403 {
404  /* TODO: lgCMatrix & igC* currently go via double prec. cholmod -
405  * since cholmod (& cs) lacks sparse 'int' matrices */
406  CHM_SP chx = AS_CHM_SP__(x);
407  int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
408  CHM_SP chxt = cholmod_transpose(chx, chx->xtype, &c);
409  SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp;
410  int tr = asLogical(tri);
411  R_CheckStack();
412 
413  tmp = VECTOR_ELT(dn, 0); /* swap the dimnames */
414  SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1));
415  SET_VECTOR_ELT(dn, 1, tmp);
416  if(!isNull(tmp = getAttrib(dn, R_NamesSymbol))) { // swap names(dimnames(.)):
417  SEXP nms_dns = PROTECT(allocVector(VECSXP, 2));
418  SET_VECTOR_ELT(nms_dns, 1, STRING_ELT(tmp, 0));
419  SET_VECTOR_ELT(nms_dns, 0, STRING_ELT(tmp, 1));
420  setAttrib(dn, R_NamesSymbol, nms_dns);
421  UNPROTECT(1);
422  }
423  UNPROTECT(1);
424  return chm_sparse_to_SEXP(chxt, 1, /* SWAP 'uplo' for triangular */
425  tr ? ((*uplo_P(x) == 'U') ? -1 : 1) : 0,
426  Rkind, tr ? diag_P(x) : "", dn);
427 }
428 
446 SEXP Csparse_Csparse_prod(SEXP a, SEXP b, SEXP bool_arith)
447 {
448  CHM_SP
449  cha = AS_CHM_SP(a),
450  chb = AS_CHM_SP(b), chc;
451  R_CheckStack();
452  static const char *valid_tri[] = { MATRIX_VALID_tri_Csparse, "" };
453  char diag[] = {'\0', '\0'};
454  int uploT = 0, nprot = 1,
455  do_bool = asLogical(bool_arith); // TRUE / NA / FALSE
456  Rboolean
457  a_is_n = (cha->xtype == CHOLMOD_PATTERN),
458  b_is_n = (chb->xtype == CHOLMOD_PATTERN),
459  force_num = (do_bool == FALSE),
460  maybe_bool= (do_bool == NA_LOGICAL);
461 
462 #ifdef DEBUG_Matrix_verbose
463  Rprintf("DBG Csparse_C*_prod(%s, %s)\n", class_P(a), class_P(b));
464 #endif
465 
466  if(a_is_n && (force_num || (maybe_bool && !b_is_n))) {
467  /* coerce 'a' to double;
468  * have no CHOLMOD function (pattern -> logical) --> use "our" code */
469  SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++;
470  cha = AS_CHM_SP(da);
471  R_CheckStack();
472  a_is_n = FALSE;
473  }
474  else if(b_is_n && (force_num || (maybe_bool && !a_is_n))) {
475  // coerce 'b' to double
476  SEXP db = PROTECT(nz2Csparse(b, x_double)); nprot++;
477  chb = AS_CHM_SP(db);
478  R_CheckStack();
479  b_is_n = FALSE;
480  }
481  chc = cholmod_ssmult(cha, chb, /*out_stype:*/ 0,
482  /* values : */ do_bool != TRUE,
483  /* sorted = TRUE: */ 1, &c);
484 
485  /* Preserve triangularity and even unit-triangularity if appropriate.
486  * Note that in that case, the multiplication itself should happen
487  * faster. But there's no support for that in CHOLMOD */
488 
489  if(R_check_class_etc(a, valid_tri) >= 0 &&
490  R_check_class_etc(b, valid_tri) >= 0)
491  if(*uplo_P(a) == *uplo_P(b)) { /* both upper, or both lower tri. */
492  uploT = (*uplo_P(a) == 'U') ? 1 : -1;
493  if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */
494  /* "remove the diagonal entries": */
495  chm_diagN2U(chc, uploT, /* do_realloc */ FALSE);
496  diag[0]= 'U';
497  }
498  else diag[0]= 'N';
499  }
500 
501  SEXP dn = PROTECT(allocVector(VECSXP, 2));
502  SET_VECTOR_ELT(dn, 0, /* establish dimnames */
503  duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0)));
504  SET_VECTOR_ELT(dn, 1,
505  duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1)));
506  UNPROTECT(nprot);
507  return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn);
508 }
509 
520 SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b, SEXP trans, SEXP bool_arith)
521 {
522  int tr = asLogical(trans), nprot = 1,
523  do_bool = asLogical(bool_arith); // TRUE / NA / FALSE
524  CHM_SP
525  cha = AS_CHM_SP(a),
526  chb = AS_CHM_SP(b),
527  chTr, chc;
528  R_CheckStack();
529  static const char *valid_tri[] = { MATRIX_VALID_tri_Csparse, "" };
530  char diag[] = {'\0', '\0'};
531  int uploT = 0;
532  Rboolean
533  a_is_n = (cha->xtype == CHOLMOD_PATTERN),
534  b_is_n = (chb->xtype == CHOLMOD_PATTERN),
535  force_num = (do_bool == FALSE),
536  maybe_bool= (do_bool == NA_LOGICAL);
537 
538  if(a_is_n && (force_num || (maybe_bool && !b_is_n))) {
539  // coerce 'a' to double
540  SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++;
541  cha = AS_CHM_SP(da);
542  R_CheckStack();
543  // a_is_n = FALSE;
544  }
545  else if(b_is_n && (force_num || (maybe_bool && !a_is_n))) {
546  // coerce 'b' to double
547  SEXP db = PROTECT(nz2Csparse(b, x_double)); nprot++;
548  chb = AS_CHM_SP(db);
549  R_CheckStack();
550  // b_is_n = FALSE;
551  }
552  else if(do_bool == TRUE) { // Want boolean arithmetic: sufficient if *one* is pattern:
553  if(!a_is_n && !b_is_n) {
554  // coerce 'a' to pattern
555  SEXP da = PROTECT(Csparse2nz(a, /* tri = */
556  R_check_class_etc(a, valid_tri) >= 0)); nprot++;
557  cha = AS_CHM_SP(da);
558  R_CheckStack();
559  // a_is_n = TRUE;
560  }
561  }
562  chTr = cholmod_transpose((tr) ? chb : cha, chb->xtype, &c);
563  chc = cholmod_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb,
564  /*out_stype:*/ 0, /* values : */ do_bool != TRUE,
565  /* sorted = TRUE: */ 1, &c);
566  cholmod_free_sparse(&chTr, &c);
567 
568  /* Preserve triangularity and unit-triangularity if appropriate;
569  * see Csparse_Csparse_prod() for comments */
570  if(R_check_class_etc(a, valid_tri) >= 0 &&
571  R_check_class_etc(b, valid_tri) >= 0)
572  if(*uplo_P(a) != *uplo_P(b)) { /* one 'U', the other 'L' */
573  uploT = (*uplo_P(b) == 'U') ? 1 : -1;
574  if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */
575  chm_diagN2U(chc, uploT, /* do_realloc */ FALSE);
576  diag[0]= 'U';
577  }
578  else diag[0]= 'N';
579  }
580 
581  SEXP dn = PROTECT(allocVector(VECSXP, 2));
582  SET_VECTOR_ELT(dn, 0, /* establish dimnames */
583  duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym),
584  (tr) ? 0 : 1)));
585  SET_VECTOR_ELT(dn, 1,
586  duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym),
587  (tr) ? 0 : 1)));
588  UNPROTECT(nprot);
589  return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn);
590 }
591 
621 SEXP Csp_dense_products(SEXP a, SEXP b,
622  Rboolean transp_a, Rboolean transp_b, Rboolean transp_ans)
623 {
624  CHM_SP cha = AS_CHM_SP(a);
625  int a_nc = transp_a ? cha->nrow : cha->ncol,
626  a_nr = transp_a ? cha->ncol : cha->nrow;
627  Rboolean
628  maybe_transp_b = (a_nc == 1),
629  b_is_vector = FALSE;
630  /* NOTE: trans_b {<--> "use t(b) instead of b" }
631  ---- "interferes" with the case automatic treatment of *vector* b.
632  In that case, t(b) or b is used "whatever make more sense",
633  according to the general R philosophy of treating vectors in matrix products.
634  */
635 
636  /* repeating a "cheap part" of mMatrix_as_dgeMatrix2(b, .) to see if
637  * we have a vector that we might 'transpose_if_vector' : */
638  static const char *valid[] = {"_NOT_A_CLASS_", MATRIX_VALID_ddense, ""};
639  /* int ctype = R_check_class_etc(b, valid);
640  * if (ctype > 0) /.* a ddenseMatrix object */
641  if (R_check_class_etc(b, valid) < 0) {
642  // not a ddenseM*: is.matrix() or vector:
643  b_is_vector = !isMatrix(b);
644  }
645 
646  if(b_is_vector) {
647  /* determine *if* we want/need to transpose at all:
648  * if (length(b) == ncol(A)) have match: use dim = c(n, 1) (<=> do *not* transp);
649  * otherwise, try to transpose: ok if (ncol(A) == 1) [see also above]: */
650  maybe_transp_b = (LENGTH(b) != a_nc);
651  // Here, we transpose already in mMatrix_as_dge*() ==> don't do it later:
652  transp_b = FALSE;
653  }
654  SEXP b_M = PROTECT(mMatrix_as_dgeMatrix2(b, maybe_transp_b));
655 
656  CHM_DN chb = AS_CHM_DN(b_M), b_t;
657  R_CheckStack();
658  int ncol_b;
659  if(transp_b) { // transpose b:
660  b_t = cholmod_allocate_dense(chb->ncol, chb->nrow, chb->ncol, chb->xtype, &c);
661  chm_transpose_dense(b_t, chb);
662  ncol_b = b_t->ncol;
663  } else
664  ncol_b = chb->ncol;
665  // Result C {with dim() before it may be transposed}:
666  CHM_DN chc = cholmod_allocate_dense(a_nr, ncol_b, a_nr, chb->xtype, &c);
667  double one[] = {1,0}, zero[] = {0,0};
668  int nprot = 2;
669 
670  /* Tim Davis, please FIXME: currently (2010-11) *fails* when a is a pattern matrix:*/
671  if(cha->xtype == CHOLMOD_PATTERN) {
672  /* warning(_("Csparse_dense_prod(): cholmod_sdmult() not yet implemented for pattern./ ngCMatrix" */
673  /* " --> slightly inefficient coercion")); */
674 
675  // This *fails* to produce a CHOLMOD_REAL ..
676  // CHM_SP chd = cholmod_l_copy(cha, cha->stype, CHOLMOD_REAL, &c);
677  // --> use our Matrix-classes
678  SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++;
679  cha = AS_CHM_SP(da);
680  }
681 
682  /* cholmod_sdmult(A, transp, alpha, beta, X, Y, &c): depending on transp == 0 / != 0:
683  * Y := alpha*(A*X) + beta*Y or alpha*(A'*X) + beta*Y; here, alpha = 1, beta = 0:
684  * Y := A*X or A'*X
685  * NB: always <sparse> %*% <dense> !
686  */
687  cholmod_sdmult(cha, transp_a, one, zero, (transp_b ? b_t : chb), /* -> */ chc, &c);
688 
689  SEXP dn = PROTECT(allocVector(VECSXP, 2)); /* establish dimnames */
690  SET_VECTOR_ELT(dn, transp_ans ? 1 : 0,
691  duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), transp_a ? 1 : 0)));
692  SET_VECTOR_ELT(dn, transp_ans ? 0 : 1,
693  duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym),
694  transp_b ? 0 : 1)));
695  if(transp_b) cholmod_free_dense(&b_t, &c);
696  UNPROTECT(nprot);
697  return chm_dense_to_SEXP(chc, 1, 0, dn, transp_ans);
698 }
699 
700 
701 SEXP Csparse_dense_prod(SEXP a, SEXP b, SEXP transp)
702 {
703  return
704  Csp_dense_products(a, b,
705  /* transp_a = */ FALSE,
706  /* transp_b = */ (*CHAR(asChar(transp)) == '2' || *CHAR(asChar(transp)) == 'B'),
707  /* transp_ans = */ (*CHAR(asChar(transp)) == 'c' || *CHAR(asChar(transp)) == 'B'));
708 }
709 
710 SEXP Csparse_dense_crossprod(SEXP a, SEXP b, SEXP transp)
711 {
712  return
713  Csp_dense_products(a, b,
714  /* transp_a = */ TRUE,
715  /* transp_b = */ (*CHAR(asChar(transp)) == '2' || *CHAR(asChar(transp)) == 'B'),
716  /* transp_ans = */ (*CHAR(asChar(transp)) == 'c' || *CHAR(asChar(transp)) == 'B'));
717 }
718 
719 
723 SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet, SEXP bool_arith)
724 {
725  int tripl = asLogical(triplet),
726  tr = asLogical(trans), /* gets reversed because _aat is tcrossprod */
727  do_bool = asLogical(bool_arith); // TRUE / NA / FALSE
728 #ifdef AS_CHM_DIAGU2N_FIXED_FINALLY
729  CHM_TR cht = tripl ? AS_CHM_TR(x) : (CHM_TR) NULL; int nprot = 1;
730 #else /* workaround needed:*/
731  SEXP xx = PROTECT(Tsparse_diagU2N(x));
732  CHM_TR cht = tripl ? AS_CHM_TR__(xx) : (CHM_TR) NULL; int nprot = 2;
733 #endif
734  CHM_SP chcp, chxt, chxc,
735  chx = (tripl ?
736  cholmod_triplet_to_sparse(cht, cht->nnz, &c) :
737  AS_CHM_SP(x));
738  SEXP dn = PROTECT(allocVector(VECSXP, 2));
739  R_CheckStack();
740  Rboolean
741  x_is_n = (chx->xtype == CHOLMOD_PATTERN),
742  x_is_sym = chx->stype != 0,
743  force_num = (do_bool == FALSE);
744 
745  if(x_is_n && force_num) {
746  // coerce 'x' to double
747  SEXP dx = PROTECT(nz2Csparse(x, x_double)); nprot++;
748  chx = AS_CHM_SP(dx);
749  R_CheckStack();
750  }
751  else if(do_bool == TRUE && !x_is_n) { // Want boolean arithmetic; need patter[n]
752  // coerce 'x' to pattern
753  static const char *valid_tri[] = { MATRIX_VALID_tri_Csparse, "" };
754  SEXP dx = PROTECT(Csparse2nz(x, /* tri = */
755  R_check_class_etc(x, valid_tri) >= 0)); nprot++;
756  chx = AS_CHM_SP(dx);
757  R_CheckStack();
758  }
759 
760  if (!tr) chxt = cholmod_transpose(chx, chx->xtype, &c);
761 
762  if (x_is_sym) // cholmod_aat() does not like symmetric
763  chxc = cholmod_copy(tr ? chx : chxt, /* stype: */ 0,
764  chx->xtype, &c);
765  // CHOLMOD/Core/cholmod_aat.c :
766  chcp = cholmod_aat(x_is_sym ? chxc : (tr ? chx : chxt),
767  (int *) NULL, 0, /* mode: */ chx->xtype, &c);
768  if(!chcp) {
769  UNPROTECT(1);
770  error(_("Csparse_crossprod(): error return from cholmod_aat()"));
771  }
772  cholmod_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c);
773  chcp->stype = 1; // symmetric
774  if (tripl) cholmod_free_sparse(&chx, &c);
775  if (!tr) cholmod_free_sparse(&chxt, &c);
776  SET_VECTOR_ELT(dn, 0, /* establish dimnames */
777  duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym),
778  (tr) ? 0 : 1)));
779  SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(dn, 0)));
780  UNPROTECT(nprot);
781  // FIXME: uploT for symmetric ?
782  return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn);
783 }
784 
787 SEXP Csparse_drop(SEXP x, SEXP tol)
788 {
789  const char *cl = class_P(x);
790  /* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */
791  int tr = (cl[1] == 't'); // FIXME - rather R_check_class_etc(..)
792  CHM_SP chx = AS_CHM_SP__(x);
793  CHM_SP ans = cholmod_copy(chx, chx->stype, chx->xtype, &c);
794  double dtol = asReal(tol);
795  int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
796  R_CheckStack();
797 
798  if(!cholmod_drop(dtol, ans, &c))
799  error(_("cholmod_drop() failed"));
800  return chm_sparse_to_SEXP(ans, 1,
801  tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
802  Rkind, tr ? diag_P(x) : "",
803  GET_SLOT(x, Matrix_DimNamesSym));
804 }
805 
808 SEXP Csparse_horzcat(SEXP x, SEXP y)
809 {
810 #define CSPARSE_CAT(_KIND_) \
811  CHM_SP chx = AS_CHM_SP__(x), chy = AS_CHM_SP__(y); \
812  R_CheckStack(); \
813  int Rk_x = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : -3, \
814  Rk_y = (chy->xtype != CHOLMOD_PATTERN) ? Real_kind(y) : -3, Rkind; \
815  if(Rk_x == -3 || Rk_y == -3) { /* at least one of them is patter"n" */ \
816  if(Rk_x == -3 && Rk_y == -3) { /* fine */ \
817  } else { /* only one is a patter"n" \
818  * "Bug" in cholmod_horzcat()/vertcat(): returns patter"n" matrix if one of them is */ \
819  Rboolean ok; \
820  if(Rk_x == -3) { \
821  ok = chm_MOD_xtype(CHOLMOD_REAL, chx, &c); Rk_x = 0; \
822  } else if(Rk_y == -3) { \
823  ok = chm_MOD_xtype(CHOLMOD_REAL, chy, &c); Rk_y = 0; \
824  } else \
825  error(_("Impossible Rk_x/Rk_y in Csparse_%s(), please report"), _KIND_); \
826  if(!ok) \
827  error(_("chm_MOD_xtype() was not successful in Csparse_%s(), please report"), \
828  _KIND_); \
829  } \
830  } \
831  Rkind = /* logical if both x and y are */ (Rk_x == 1 && Rk_y == 1) ? 1 : 0
832 
833  CSPARSE_CAT("horzcat");
834  // TODO: currently drops dimnames - and we fix at R level;
835 
836  return chm_sparse_to_SEXP(cholmod_horzcat(chx, chy, 1, &c),
837  1, 0, Rkind, "", R_NilValue);
838 }
839 
842 SEXP Csparse_vertcat(SEXP x, SEXP y)
843 {
844  CSPARSE_CAT("vertcat");
845  // TODO: currently drops dimnames - and we fix at R level;
846 
847  return chm_sparse_to_SEXP(cholmod_vertcat(chx, chy, 1, &c),
848  1, 0, Rkind, "", R_NilValue);
849 }
850 
851 SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2)
852 {
853  CHM_SP chx = AS_CHM_SP__(x);
854  int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
855  CHM_SP ans = cholmod_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c);
856  R_CheckStack();
857 
858  return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "",
859  GET_SLOT(x, Matrix_DimNamesSym));
860 }
861 
862 SEXP Csparse_diagU2N(SEXP x)
863 {
864  const char *cl = class_P(x);
865  /* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */
866  if (cl[1] != 't' || *diag_P(x) != 'U') {
867  /* "trivially fast" when not triangular (<==> no 'diag' slot),
868  or not *unit* triangular */
869  return (x);
870  }
871  else { /* unit triangular (diag='U'): "fill the diagonal" & diag:= "N" */
872  CHM_SP chx = AS_CHM_SP__(x);
873  CHM_SP eye = cholmod_speye(chx->nrow, chx->ncol, chx->xtype, &c);
874  double one[] = {1, 0};
875  CHM_SP ans = cholmod_add(chx, eye, one, one, TRUE, TRUE, &c);
876  int uploT = (*uplo_P(x) == 'U') ? 1 : -1;
877  int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
878 
879  R_CheckStack();
880  cholmod_free_sparse(&eye, &c);
881  return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N",
882  GET_SLOT(x, Matrix_DimNamesSym));
883  }
884 }
885 
886 SEXP Csparse_diagN2U(SEXP x)
887 {
888  const char *cl = class_P(x);
889  /* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */
890  if (cl[1] != 't' || *diag_P(x) != 'N') {
891  /* "trivially fast" when not triangular (<==> no 'diag' slot),
892  or already *unit* triangular */
893  return (x);
894  }
895  else { /* triangular with diag='N'): now drop the diagonal */
896  /* duplicate, since chx will be modified: */
897  SEXP xx = PROTECT(duplicate(x));
898  CHM_SP chx = AS_CHM_SP__(xx);
899  int uploT = (*uplo_P(x) == 'U') ? 1 : -1,
900  Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
901  R_CheckStack();
902 
903  chm_diagN2U(chx, uploT, /* do_realloc */ FALSE);
904 
905  SEXP ans = chm_sparse_to_SEXP(chx, /*dofree*/ 0/* or 1 ?? */,
906  uploT, Rkind, "U",
907  GET_SLOT(x, Matrix_DimNamesSym));
908  UNPROTECT(1);// only now !
909  return ans;
910  }
911 }
912 
922 SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j)
923 {
924  CHM_SP chx = AS_CHM_SP(x); /* << does diagU2N() when needed */
925  int rsize = (isNull(i)) ? -1 : LENGTH(i),
926  csize = (isNull(j)) ? -1 : LENGTH(j);
927  int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
928  R_CheckStack();
929 
930  if (rsize >= 0 && !isInteger(i))
931  error(_("Index i must be NULL or integer"));
932  if (csize >= 0 && !isInteger(j))
933  error(_("Index j must be NULL or integer"));
934 
935 #define CHM_SUB(_M_, _i_, _j_) \
936  cholmod_submatrix(_M_, \
937  (rsize < 0) ? NULL : INTEGER(_i_), rsize, \
938  (csize < 0) ? NULL : INTEGER(_j_), csize, \
939  TRUE, TRUE, &c)
940  CHM_SP ans;
941  if (!chx->stype) {/* non-symmetric Matrix */
942  ans = CHM_SUB(chx, i, j);
943  }
944  else {
945  /* for now, cholmod_submatrix() only accepts "generalMatrix" */
946  CHM_SP tmp = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);
947  ans = CHM_SUB(tmp, i, j);
948  cholmod_free_sparse(&tmp, &c);
949  }
950 
951  // "FIXME": currently dropping dimnames, and adding them afterwards in R :
952  /* // dimnames: */
953  /* SEXP x_dns = GET_SLOT(x, Matrix_DimNamesSym), */
954  /* dn = PROTECT(allocVector(VECSXP, 2)); */
955  return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", /* dimnames: */ R_NilValue);
956 }
957 
958 #define _d_Csp_
959 #include "t_Csparse_subassign.c"
960 
961 #define _l_Csp_
962 #include "t_Csparse_subassign.c"
963 
964 #define _i_Csp_
965 #include "t_Csparse_subassign.c"
966 
967 #define _n_Csp_
968 #include "t_Csparse_subassign.c"
969 
970 #define _z_Csp_
971 #include "t_Csparse_subassign.c"
972 
973 
974 
975 SEXP Csparse_MatrixMarket(SEXP x, SEXP fname)
976 {
977  FILE *f = fopen(CHAR(asChar(fname)), "w");
978 
979  if (!f)
980  error(_("failure to open file \"%s\" for writing"),
981  CHAR(asChar(fname)));
982  if (!cholmod_write_sparse(f, AS_CHM_SP(x),
983  (CHM_SP)NULL, (char*) NULL, &c))
984  error(_("cholmod_write_sparse returned error code"));
985  fclose(f);
986  return R_NilValue;
987 }
988 
989 
1002 SEXP diag_tC_ptr(int n, int *x_p, double *x_x, Rboolean is_U, int *perm,
1003 /* ^^^^^^ FIXME[Generalize] to int / ... */
1004  SEXP resultKind)
1005 {
1006  const char* res_ch = CHAR(STRING_ELT(resultKind,0));
1007  enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log, min, max, range
1008  } res_kind = ((!strcmp(res_ch, "trace")) ? trace :
1009  ((!strcmp(res_ch, "sumLog")) ? sum_log :
1010  ((!strcmp(res_ch, "prod")) ? prod :
1011  ((!strcmp(res_ch, "min")) ? min :
1012  ((!strcmp(res_ch, "max")) ? max :
1013  ((!strcmp(res_ch, "range")) ? range :
1014  ((!strcmp(res_ch, "diag")) ? diag :
1015  ((!strcmp(res_ch, "diagBack")) ? diag_backpermuted :
1016  -1))))))));
1017  int i, n_x, i_from;
1018  SEXP ans = PROTECT(allocVector(REALSXP,
1019 /* ^^^^ FIXME[Generalize] */
1020  (res_kind == diag ||
1021  res_kind == diag_backpermuted) ? n :
1022  (res_kind == range ? 2 : 1)));
1023  double *v = REAL(ans);
1024 /* ^^^^^^ ^^^^ FIXME[Generalize] */
1025 
1026  i_from = (is_U ? -1 : 0);
1027 
1028 #define for_DIAG(v_ASSIGN) \
1029  for(i = 0; i < n; i++) { \
1030  /* looking at i-th column */ \
1031  n_x = x_p[i+1] - x_p[i];/* #{entries} in this column */ \
1032  if( is_U) i_from += n_x; \
1033  v_ASSIGN; \
1034  if(!is_U) i_from += n_x; \
1035  }
1036 
1037  /* NOTA BENE: we assume -- uplo = "L" i.e. lower triangular matrix
1038  * for uplo = "U" (makes sense with a "dtCMatrix" !),
1039  * should use x_x[i_from + (n_x - 1)] instead of x_x[i_from],
1040  * where n_x = (x_p[i+1] - x_p[i])
1041  */
1042 
1043  switch(res_kind) {
1044  case trace: // = sum
1045  v[0] = 0.;
1046  for_DIAG(v[0] += x_x[i_from]);
1047  break;
1048 
1049  case sum_log:
1050  v[0] = 0.;
1051  for_DIAG(v[0] += log(x_x[i_from]));
1052  break;
1053 
1054  case prod:
1055  v[0] = 1.;
1056  for_DIAG(v[0] *= x_x[i_from]);
1057  break;
1058 
1059  case min:
1060  v[0] = R_PosInf;
1061  for_DIAG(if(v[0] > x_x[i_from]) v[0] = x_x[i_from]);
1062  break;
1063 
1064  case max:
1065  v[0] = R_NegInf;
1066  for_DIAG(if(v[0] < x_x[i_from]) v[0] = x_x[i_from]);
1067  break;
1068 
1069  case range:
1070  v[0] = R_PosInf;
1071  v[1] = R_NegInf;
1072  for_DIAG(if(v[0] > x_x[i_from]) v[0] = x_x[i_from];
1073  if(v[1] < x_x[i_from]) v[1] = x_x[i_from]);
1074  break;
1075 
1076  case diag:
1077  for_DIAG(v[i] = x_x[i_from]);
1078  break;
1079 
1080  case diag_backpermuted:
1081  for_DIAG(v[i] = x_x[i_from]);
1082 
1083  warning(_("%s = '%s' (back-permuted) is experimental"),
1084  "resultKind", "diagBack");
1085  /* now back_permute : */
1086  for(i = 0; i < n; i++) {
1087  double tmp = v[i]; v[i] = v[perm[i]]; v[perm[i]] = tmp;
1088  /*^^^^ FIXME[Generalize] */
1089  }
1090  break;
1091 
1092  default: /* -1 from above */
1093  error(_("diag_tC(): invalid 'resultKind'"));
1094  /* Wall: */ ans = R_NilValue; v = REAL(ans);
1095  }
1096 
1097  UNPROTECT(1);
1098  return ans;
1099 }
1100 
1114 SEXP diag_tC(SEXP obj, SEXP resultKind)
1115 {
1116 
1117  SEXP
1118  pslot = GET_SLOT(obj, Matrix_pSym),
1119  xslot = GET_SLOT(obj, Matrix_xSym);
1120  Rboolean is_U = (R_has_slot(obj, Matrix_uploSym) &&
1121  *CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))) == 'U');
1122  int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */
1123  *x_p = INTEGER(pslot), pp = -1, *perm;
1124  double *x_x = REAL(xslot);
1125 /* ^^^^^^ ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/
1126 
1127  if(R_has_slot(obj, Matrix_permSym))
1128  perm = INTEGER(GET_SLOT(obj, Matrix_permSym));
1129  else perm = &pp;
1130 
1131  return diag_tC_ptr(n, x_p, x_x, is_U, perm, resultKind);
1132 }
1133 
1134 
1153 SEXP create_Csparse(char* cls, int* i, int* j, int* p, int np,
1154  void* x, int nnz, int* dims, SEXP dimnames,
1155  int index1)
1156 {
1157  SEXP ans;
1158  int *ij = (int*)NULL, *tri, *trj,
1159  mi, mj, mp, nrow = -1, ncol = -1;
1160  int xtype = -1; /* -Wall */
1161  CHM_TR T;
1162  CHM_SP A;
1163 
1164  if (np < 0 || nnz < 0)
1165  error(_("negative vector lengths not allowed: np = %d, nnz = %d"),
1166  np, nnz);
1167  if (1 != ((mi = (i == (int*)NULL)) +
1168  (mj = (j == (int*)NULL)) +
1169  (mp = (p == (int*)NULL))))
1170  error(_("exactly 1 of 'i', 'j' or 'p' must be NULL"));
1171  if (mp) {
1172  if (np) error(_("np = %d, must be zero when p is NULL"), np);
1173  } else {
1174  if (np) { /* Expand p to form i or j */
1175  if (!(p[0])) error(_("p[0] = %d, should be zero"), p[0]);
1176  for (int ii = 0; ii < np; ii++)
1177  if (p[ii] > p[ii + 1])
1178  error(_("p must be non-decreasing"));
1179  if (p[np] != nnz)
1180  error("p[np] = %d != nnz = %d", p[np], nnz);
1181  ij = Calloc(nnz, int);
1182  if (mi) {
1183  i = ij;
1184  nrow = np;
1185  } else {
1186  j = ij;
1187  ncol = np;
1188  }
1189  /* Expand p to 0-based indices */
1190  for (int ii = 0; ii < np; ii++)
1191  for (int jj = p[ii]; jj < p[ii + 1]; jj++) ij[jj] = ii;
1192  } else {
1193  if (nnz)
1194  error(_("Inconsistent dimensions: np = 0 and nnz = %d"),
1195  nnz);
1196  }
1197  }
1198  /* calculate nrow and ncol */
1199  if (nrow < 0) {
1200  for (int ii = 0; ii < nnz; ii++) {
1201  int i1 = i[ii] + (index1 ? 0 : 1); /* 1-based index */
1202  if (i1 < 1) error(_("invalid row index at position %d"), ii);
1203  if (i1 > nrow) nrow = i1;
1204  }
1205  }
1206  if (ncol < 0) {
1207  for (int jj = 0; jj < nnz; jj++) {
1208  int j1 = j[jj] + (index1 ? 0 : 1);
1209  if (j1 < 1) error(_("invalid column index at position %d"), jj);
1210  if (j1 > ncol) ncol = j1;
1211  }
1212  }
1213  if (dims != (int*)NULL) {
1214  if (dims[0] > nrow) nrow = dims[0];
1215  if (dims[1] > ncol) ncol = dims[1];
1216  }
1217  /* check the class name */
1218  if (strlen(cls) != 8)
1219  error(_("strlen of cls argument = %d, should be 8"), strlen(cls));
1220  if (!strcmp(cls + 2, "CMatrix"))
1221  error(_("cls = \"%s\" does not end in \"CMatrix\""), cls);
1222  switch(cls[0]) {
1223  case 'd':
1224  case 'l':
1225  xtype = CHOLMOD_REAL;
1226  break;
1227  case 'n':
1228  xtype = CHOLMOD_PATTERN;
1229  break;
1230  default:
1231  error(_("cls = \"%s\" must begin with 'd', 'l' or 'n'"), cls);
1232  }
1233  if (cls[1] != 'g')
1234  error(_("Only 'g'eneral sparse matrix types allowed"));
1235  /* allocate and populate the triplet */
1236  T = cholmod_allocate_triplet((size_t)nrow, (size_t)ncol, (size_t)nnz, 0,
1237  xtype, &c);
1238  T->x = x;
1239  tri = (int*)T->i;
1240  trj = (int*)T->j;
1241  for (int ii = 0; ii < nnz; ii++) {
1242  tri[ii] = i[ii] - ((!mi && index1) ? 1 : 0);
1243  trj[ii] = j[ii] - ((!mj && index1) ? 1 : 0);
1244  }
1245  /* create the cholmod_sparse structure */
1246  A = cholmod_triplet_to_sparse(T, nnz, &c);
1247  cholmod_free_triplet(&T, &c);
1248  /* copy the information to the SEXP */
1249  ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cls)));
1250 // FIXME: This has been copied from chm_sparse_to_SEXP in chm_common.c
1251  /* allocate and copy common slots */
1252  nnz = cholmod_nnz(A, &c);
1253  dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
1254  dims[0] = A->nrow; dims[1] = A->ncol;
1255  Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, A->ncol + 1)), (int*)A->p, A->ncol + 1);
1256  Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)), (int*)A->i, nnz);
1257  switch(cls[1]) {
1258  case 'd':
1259  Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)), (double*)A->x, nnz);
1260  break;
1261  case 'l':
1262  error(_("code not yet written for cls = \"lgCMatrix\""));
1263  }
1264 /* FIXME: dimnames are *NOT* put there yet (if non-NULL) */
1265  cholmod_free_sparse(&A, &c);
1266  UNPROTECT(1);
1267  return ans;
1268 }
SEXP Matrix_DimSym
Definition: Syms.h:2
#define class_P(_x_)
Definition: Mutils.h:178
SEXP chm_dense_to_SEXP(CHM_DN a, int dofree, int Rkind, SEXP dn, Rboolean transp)
Copy the contents of a to an appropriate denseMatrix object and, optionally, free a or free both a an...
Definition: chm_common.c:835
#define AS_CHM_SP__(x)
Definition: chm_common.h:49
x_slot_kind
Definition: Mutils.h:184
SEXP Csparse_drop(SEXP x, SEXP tol)
cholmod_sparse * CHM_SP
Definition: chm_common.h:25
SEXP diag_tC_ptr(int n, int *x_p, double *x_x, Rboolean is_U, int *perm, SEXP resultKind)
SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri)
Definition: Csparse.c:210
static R_INLINE SEXP mMatrix_as_dgeMatrix2(SEXP A, Rboolean tr_if_vec)
Definition: Mutils.h:338
#define MATRIX_VALID_tri_Csparse
Definition: Mutils.h:402
SEXP Matrix_xSym
Definition: Syms.h:2
#define slot_dup(dest, src, sym)
Definition: Mutils.h:149
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
cholmod_triplet * CHM_TR
Definition: chm_common.h:27
SEXP Csparse2nz(SEXP x, Rboolean tri)
Definition: Csparse.c:198
SEXP chm_sparse_to_SEXP(CHM_SP a, int dofree, int uploT, int Rkind, const char *diag, SEXP dn)
Copy the contents of a to an appropriate CsparseMatrix object and, optionally, free a or free both a ...
Definition: chm_common.c:335
SEXP Csparse_validate(SEXP x)
Definition: Csparse.c:39
cholmod_common cl
Definition: chm_common.c:16
SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b, SEXP trans, SEXP bool_arith)
[t]crossprod (, )
Definition: Csparse.c:520
SEXP Matrix_uploSym
Definition: Syms.h:2
SEXP Csparse_to_Tsparse(SEXP x, SEXP tri)
Definition: Csparse.c:305
SEXP Csp_dense_products(SEXP a, SEXP b, Rboolean transp_a, Rboolean transp_b, Rboolean transp_ans)
SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j)
SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet, SEXP bool_arith)
SEXP Csparse_to_matrix(SEXP x, SEXP chk, SEXP symm)
Definition: Csparse.c:284
SEXP Csparse_MatrixMarket(SEXP x, SEXP fname)
#define AS_CHM_TR__(x)
Definition: chm_common.h:50
SEXP nz2Csparse(SEXP x, enum x_slot_kind r_kind)
Definition: Csparse.c:228
Rboolean equal_string_vectors(SEXP s1, SEXP s2)
Definition: Mutils.c:286
SEXP Csparse_to_vector(SEXP x)
Definition: Csparse.c:300
SEXP Csparse_validate_(SEXP x, Rboolean maybe_modify)
SEXP Matrix_jSym
Definition: Syms.h:2
SEXP Csparse_sort(SEXP x)
Definition: Csparse.c:56
static int xtype(int ctype)
xtype: the kind of numeric (think "x slot") of Cholmod sparse matrices.
Definition: chm_common.c:126
SEXP Csparse_horzcat(SEXP x, SEXP y)
SEXP Csparse_symmetric_to_general(SEXP x)
Definition: Csparse.c:343
SEXP chm_dense_to_vector(CHM_DN a, int dofree)
Copy the contents of a to a numeric R object and, optionally, free a or free both a and its pointer t...
Definition: chm_common.c:996
SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo, SEXP sym_dmns)
Definition: Csparse.c:357
SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2)
#define MATRIX_VALID_ddense
Definition: Mutils.h:360
SEXP Matrix_permSym
Definition: Syms.h:2
#define uplo_P(_x_)
Definition: Mutils.h:174
SEXP Csparse_Csparse_prod(SEXP a, SEXP b, SEXP bool_arith)
A %*% B - for matrices of class CsparseMatrix (R package "Matrix")
Definition: Csparse.c:446
#define _(String)
Definition: Mutils.h:32
cholmod_dense * CHM_DN
Definition: chm_common.h:21
SEXP Matrix_iSym
Definition: Syms.h:2
SEXP chm_dense_to_matrix(CHM_DN a, int dofree, SEXP dn)
Copy the contents of a to a matrix object and, optionally, free a or free both a and its pointer to i...
Definition: chm_common.c:942
SEXP Csparse_to_tTsparse(SEXP x, SEXP uplo, SEXP diag)
Definition: Csparse.c:330
SEXP create_Csparse(char *cls, int *i, int *j, int *p, int np, void *x, int nnz, int *dims, SEXP dimnames, int index1)
#define AS_CHM_DN(x)
Definition: chm_common.h:43
#define slot_dup_if_has(dest, src, sym)
Definition: Mutils.h:156
#define diag_P(_x_)
Definition: Mutils.h:175
SEXP diag_tC(SEXP obj, SEXP resultKind)
#define Real_kind(_x_)
Definition: Mutils.h:195
SEXP Matrix_pSym
Definition: Syms.h:2
SEXP Csparse_validate2(SEXP x, SEXP maybe_modify)
Definition: Csparse.c:51
#define AS_CHM_SP2(x, chk)
Definition: chm_common.h:52
SEXP Csparse_diagN2U(SEXP x)
SEXP symmetric_DimNames(SEXP dn)
Produce symmetric 'Dimnames' from possibly asymmetric ones.
Definition: Mutils.c:1269
SEXP Tsparse_diagU2N(SEXP x)
Definition: Tsparse.c:63
SEXP Csparse_transpose(SEXP x, SEXP tri)
Definition: Csparse.c:402
SEXP Matrix_diagSym
Definition: Syms.h:2
void chm_diagN2U(CHM_SP chx, int uploT, Rboolean do_realloc)
Drop the (unit) diagonal entries from a cholmod_sparse matrix.
Definition: chm_common.c:1191
#define MATRIX_VALID_nCsparse
Definition: Mutils.h:381
SEXP nz_pattern_to_Csparse(SEXP x, SEXP res_kind)
Definition: Csparse.c:221
SEXP Csparse_diagU2N(SEXP x)
SEXP Csparse_to_dense(SEXP x, SEXP symm_or_tri)
From a CsparseMatrix, produce a dense one.
Definition: Csparse.c:121
SEXP Csparse_to_tCsparse(SEXP x, SEXP uplo, SEXP diag)
Definition: Csparse.c:319
cholmod_common c
Definition: chm_common.c:15
Rboolean isValid_Csparse(SEXP x)
"Cheap" C version of Csparse_validate() - not sorting :
Definition: Csparse.c:11
#define MATRIX_VALID_Csparse
Definition: Mutils.h:384
SEXP Rsparse_validate(SEXP x)
Definition: Csparse.c:62
SEXP Csparse_vertcat(SEXP x, SEXP y)
SEXP Csparse_dense_prod(SEXP a, SEXP b, SEXP transp)
static int is_sym(cs *A)
Definition: cs_utils.c:5
SEXP chm_triplet_to_SEXP(CHM_TR a, int dofree, int uploT, int Rkind, const char *diag, SEXP dn)
Copy the contents of a to an appropriate TsparseMatrix object and, optionally, free a or free both a ...
Definition: chm_common.c:573
SEXP Csparse_dense_crossprod(SEXP a, SEXP b, SEXP transp)
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
#define AS_CHM_SP(x)
Definition: chm_common.h:46
#define AS_CHM_TR(x)
Definition: chm_common.h:47
void chm_transpose_dense(CHM_DN ans, CHM_DN x)
Transpose a cholmod_dense matrix ("too trivial to be in CHOLMOD?")
Definition: chm_common.c:754