Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
Tsparse.c
Go to the documentation of this file.
1  /* Sparse matrices in triplet form */
2 #include "Tsparse.h"
3 #include "chm_common.h"
4 
5 SEXP Tsparse_validate(SEXP x)
6 {
7  /* NB: we do *NOT* check a potential 'x' slot here, at all */
8  SEXP
9  islot = GET_SLOT(x, Matrix_iSym),
10  jslot = GET_SLOT(x, Matrix_jSym),
11  dimslot = GET_SLOT(x, Matrix_DimSym);
12  int j,
13  nrow = INTEGER(dimslot)[0],
14  ncol = INTEGER(dimslot)[1],
15  nnz = length(islot),
16  *xj = INTEGER(jslot),
17  *xi = INTEGER(islot);
18 
19  if (length(jslot) != nnz)
20  return mkString(_("lengths of slots i and j must match"));
21  /* FIXME: this is checked in super class -- no need to do here: */
22  if (length(dimslot) != 2)
23  return mkString(_("slot Dim must have length 2"));
24 
25  for (j = 0; j < nnz; j++) {
26  if (xi[j] < 0 || xi[j] >= nrow)
27  return mkString(_("all row indices (slot 'i') must be between 0 and nrow-1 in a TsparseMatrix"));
28  if (xj[j] < 0 || xj[j] >= ncol)
29  return mkString(_("all column indices (slot 'j') must be between 0 and ncol-1 in a TsparseMatrix"));
30  }
31  return ScalarLogical(1);
32 }
33 
34 SEXP Tsparse_to_Csparse(SEXP x, SEXP tri)
35 {
36  CHM_TR chxt = AS_CHM_TR__(x); /* << should *preserve* diag = "U" ! */
37  CHM_SP chxs = cholmod_triplet_to_sparse(chxt, chxt->nnz, &c);
38  int tr = asLogical(tri);
39  int Rkind = (chxt->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
40  R_CheckStack();
41 
42  return chm_sparse_to_SEXP(chxs, 1,
43  tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
44  Rkind, tr ? diag_P(x) : "",
45  GET_SLOT(x, Matrix_DimNamesSym));
46 }
47 
48 /* speedup utility, needed e.g. after subsetting: */
49 SEXP Tsparse_to_tCsparse(SEXP x, SEXP uplo, SEXP diag)
50 {
51  CHM_TR chxt = AS_CHM_TR__(x);
52  CHM_SP chxs = cholmod_triplet_to_sparse(chxt, chxt->nnz, &c);
53  int Rkind = (chxt->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
54  R_CheckStack();
55 
56  return chm_sparse_to_SEXP(chxs, 1,
57  /* uploT = */ (*CHAR(asChar(uplo)) == 'U')? 1: -1,
58  Rkind,
59  /* diag = */ CHAR(STRING_ELT(diag, 0)),
60  GET_SLOT(x, Matrix_DimNamesSym));
61 }
62 
63 SEXP Tsparse_diagU2N(SEXP x)
64 {
65  static const char *valid[] = {
66  "dtTMatrix", /* 0 */
67  "ltTMatrix", /* 1 */
68  "ntTMatrix", /* 2 : no x slot */
69  "ztTMatrix", /* 3 */
70  ""};
71 /* #define xSXP(iTyp) ((iTyp == 0) ? REALSXP : ((iTyp == 1) ? LGLSXP : /\* else *\/ CPLXSXP)); */
72 /* #define xTYPE(iTyp) ((iTyp == 0) ? double : ((iTyp == 1) ? int : /\* else *\/ Rcomplex)); */
73  int ctype = R_check_class_etc(x, valid);
74 
75  if (ctype < 0 || *diag_P(x) != 'U') {
76  /* "trivially fast" when not triangular (<==> no 'diag' slot),
77  or not *unit* triangular */
78  return (x);
79  }
80  else { /* instead of going to Csparse -> Cholmod -> Csparse -> Tsparse, work directly: */
81  int i, n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0],
82  nnz = length(GET_SLOT(x, Matrix_iSym)),
83  new_n = nnz + n;
84  SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS(class_P(x))));
85  int *islot = INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, new_n)),
86  *jslot = INTEGER(ALLOC_SLOT(ans, Matrix_jSym, INTSXP, new_n));
87 
88  slot_dup(ans, x, Matrix_DimSym);
89  SET_DimNames(ans, x);
90  slot_dup(ans, x, Matrix_uploSym);
91  SET_SLOT(ans, Matrix_diagSym, mkString("N"));
92 
93  /* Build the new i- and j- slots : first copy the current : */
94  Memcpy(islot, INTEGER(GET_SLOT(x, Matrix_iSym)), nnz);
95  Memcpy(jslot, INTEGER(GET_SLOT(x, Matrix_jSym)), nnz);
96  /* then, add the new (i,j) slot entries: */
97  for(i = 0; i < n; i++) {
98  islot[i + nnz] = i;
99  jslot[i + nnz] = i;
100  }
101 
102  /* build the new x-slot : */
103  switch(ctype) {
104  case 0: { /* "d" */
105  double *x_new = REAL(ALLOC_SLOT(ans, Matrix_xSym,
106  REALSXP, new_n));
107  Memcpy(x_new, REAL(GET_SLOT(x, Matrix_xSym)), nnz);
108  for(i = 0; i < n; i++) /* add x[i,i] = 1. */
109  x_new[i + nnz] = 1.;
110  break;
111  }
112  case 1: { /* "l" */
113  int *x_new = LOGICAL(ALLOC_SLOT(ans, Matrix_xSym,
114  LGLSXP, new_n));
115  Memcpy(x_new, LOGICAL(GET_SLOT(x, Matrix_xSym)), nnz);
116  for(i = 0; i < n; i++) /* add x[i,i] = 1 (= TRUE) */
117  x_new[i + nnz] = 1;
118  break;
119  }
120  case 2: /* "n" */
121  /* nothing to do here */
122  break;
123 
124  case 3: { /* "z" */
125  Rcomplex *x_new = COMPLEX(ALLOC_SLOT(ans, Matrix_xSym,
126  CPLXSXP, new_n));
127  Memcpy(x_new, COMPLEX(GET_SLOT(x, Matrix_xSym)), nnz);
128  for(i = 0; i < n; i++) /* add x[i,i] = 1 (= TRUE) */
129  x_new[i + nnz] = (Rcomplex) {1., 0.};
130  break;
131  }
132 
133  }/* switch() */
134 
135  UNPROTECT(1);
136  return ans;
137  }
138 }
SEXP Matrix_DimSym
Definition: Syms.h:2
#define class_P(_x_)
Definition: Mutils.h:178
cholmod_sparse * CHM_SP
Definition: chm_common.h:25
SEXP Matrix_xSym
Definition: Syms.h:2
#define slot_dup(dest, src, sym)
Definition: Mutils.h:149
SEXP Tsparse_to_Csparse(SEXP x, SEXP tri)
Definition: Tsparse.c:34
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
cholmod_triplet * CHM_TR
Definition: chm_common.h:27
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 Matrix_uploSym
Definition: Syms.h:2
#define AS_CHM_TR__(x)
Definition: chm_common.h:50
SEXP Matrix_jSym
Definition: Syms.h:2
#define uplo_P(_x_)
Definition: Mutils.h:174
#define _(String)
Definition: Mutils.h:32
SEXP Matrix_iSym
Definition: Syms.h:2
#define diag_P(_x_)
Definition: Mutils.h:175
#define Real_kind(_x_)
Definition: Mutils.h:195
SEXP Tsparse_diagU2N(SEXP x)
Definition: Tsparse.c:63
SEXP Matrix_diagSym
Definition: Syms.h:2
cholmod_common c
Definition: chm_common.c:15
static R_INLINE void SET_DimNames(SEXP dest, SEXP src)
Definition: Mutils.h:161
SEXP Tsparse_validate(SEXP x)
Definition: Tsparse.c:5
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
SEXP Tsparse_to_tCsparse(SEXP x, SEXP uplo, SEXP diag)
Definition: Tsparse.c:49