Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
cs_utils.c
Go to the documentation of this file.
1 #include "cs_utils.h"
2 
3 /* Borrowed from one of Tim Davis' examples in the CSparse Demo directory */
4 /* 1 if A is square & upper tri., -1 if square & lower tri., 0 otherwise */
5 static int is_sym (cs *A)
6 {
7  int is_upper, is_lower, j, p, n = A->n, m = A->m, *Ap = A->p, *Ai = A->i ;
8  if (m != n) return (0) ;
9  is_upper = 1 ;
10  is_lower = 1 ;
11  for (j = 0 ; j < n ; j++)
12  {
13  for (p = Ap [j] ; p < Ap [j+1] ; p++)
14  {
15  if (Ai [p] > j) is_upper = 0 ;
16  if (Ai [p] < j) is_lower = 0 ;
17  }
18  }
19  return (is_upper ? 1 : (is_lower ? -1 : 0)) ;
20 }
21 
22 
31 static CSP csp_eye(int n)
32 {
33  CSP eye = cs_spalloc(n, n, n, 1, 0);
34  int *ep = eye->p, *ei = eye->i;
35  double *ex = eye->x;
36 
37  if (n <= 0) error(_("csp_eye argument n must be positive"));
38  eye->nz = -1; /* compressed column storage */
39  for (int j = 0; j < n; j++) {
40  ep[j] = ei[j] = j;
41  ex[j] = 1;
42  }
43  eye->nzmax = ep[n] = n;
44  return eye;
45 }
46 
61 cs *Matrix_as_cs(cs *ans, SEXP x, Rboolean check_Udiag)
62 {
63  static const char *valid[] = {"dgCMatrix", "dtCMatrix", ""};
64  /* had also "dsCMatrix", but that only stores one triangle */
65  int *dims, ctype = R_check_class_etc(x, valid);
66  SEXP islot;
67 
68  if (ctype < 0) error(_("invalid class of 'x' in Matrix_as_cs(a, x)"));
69  /* dimensions and nzmax */
70  dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
71  ans->m = dims[0]; ans->n = dims[1];
72  islot = GET_SLOT(x, Matrix_iSym);
73  ans->nz = -1; /* indicates compressed column storage */
74  ans->nzmax = LENGTH(islot);
75  ans->i = INTEGER(islot);
76  ans->p = INTEGER(GET_SLOT(x, Matrix_pSym));
77  ans->x = REAL(GET_SLOT(x, Matrix_xSym));
78 
79  if(check_Udiag && ctype == 1 && (*diag_P(x) == 'U')) { /* diagU2N(.) : */
80  int n = dims[0];
81  CSP I_n = csp_eye(n);
82 
83  /* tmp := 1*ans + 1*eye -- result is newly allocated in cs_add(): */
84  CSP tmp = cs_add(ans, I_n, 1., 1.), t2;
85  int nz = (tmp->p)[n];
86 
87  /* double transpose trick to sort the columns */
88  cs_spfree(I_n);
89  t2 = cs_transpose(tmp, 1); /* transpose including values */
90  cs_spfree(tmp);
91  tmp = cs_transpose(t2, 1);
92  cs_spfree(t2);
93 
94  /* content(ans) := content(tmp) : */
95  ans->nzmax = nz;
96  /* The ans "slots" were pointers to x@ <slots>; all need new content now: */
97  ans->p = Memcpy(( int*) R_alloc(n+1, sizeof(int)),
98  ( int*) tmp->p, n+1);
99  ans->i = Memcpy(( int*) R_alloc(nz, sizeof(int)),
100  ( int*) tmp->i, nz);
101  ans->x = Memcpy((double*) R_alloc(nz, sizeof(double)),
102  (double*) tmp->x, nz);
103 
104  cs_spfree(tmp);
105  }
106  return ans;
107 }
108 
120 // FIXME: Change API : Use object, not just class name 'cl' -- and use R_check_class(obj, *)
121 SEXP Matrix_cs_to_SEXP(cs *a, char *cl, int dofree, SEXP dn)
122 {
123  static const char *valid[] = {"dgCMatrix", "dsCMatrix", "dtCMatrix", ""};
124  int ctype = Matrix_check_class(cl, valid);
125 
126  if (ctype < 0)
127  error(_("invalid class of object to %s"), "Matrix_cs_to_SEXP");
128  SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cl)));
129  /* allocate and copy common slots */
130  int *dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
131  PROTECT(dn); // <- as in chm_sparse_to_SEXP()
132  dims[0] = a->m; dims[1] = a->n;
133  Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, a->n + 1)),
134  a->p, a->n + 1);
135  int nz = a->p[a->n];
136  Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), a->i, nz);
137  Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), a->x, nz);
138  if (ctype > 0) { /* dsC or dtC */
139  int uplo = is_sym(a);
140  if (!uplo)
141  error(_("cs matrix not compatible with class '%s'"), valid[ctype]);
142  if (ctype == 2) /* dtC* */
143  SET_SLOT(ans, Matrix_diagSym, mkString("N"));
144  SET_SLOT(ans, Matrix_uploSym, mkString(uplo < 0 ? "L" : "U"));
145  }
146  if (dofree > 0) cs_spfree(a);
147  if (dofree < 0) Free(a);
148  if (dn != R_NilValue)
149  SET_SLOT(ans, Matrix_DimNamesSym, duplicate(dn));
150  UNPROTECT(2);
151  return ans;
152 }
153 
154 #if 0 /* unused ------------------------------------*/
155 /* -------------------------------------*/
156 
166 css *Matrix_as_css(css *ans, SEXP x)
167 {
168  char *cl = class_P(x);
169  static const char *valid[] = {"css_LU", "css_QR", ""};
170  int *nz = INTEGER(GET_SLOT(x, install("nz"))),
171  ctype = Matrix_check_class(cl, valid);
172 
173  if (ctype < 0)
174  error(_("invalid class of object to %s"), "Matrix_as_css");
175  ans->q = INTEGER(GET_SLOT(x, install("Q")));
176  ans->m2 = nz[0]; ans->lnz = nz[1]; ans->unz = nz[2];
177  switch(ctype) {
178  case 0: /* css_LU */
179  ans->pinv = (int *) NULL;
180  ans->parent = (int *) NULL;
181  ans->cp = (int *) NULL;
182  break;
183  case 1: /* css_QR */
184  ans->pinv = INTEGER(GET_SLOT(x, install("Pinv")));
185  ans->parent = INTEGER(GET_SLOT(x, install("parent")));
186  ans->cp = INTEGER(GET_SLOT(x, install("cp")));
187  break;
188  default:
189  error(_("invalid class of object to %s"), "Matrix_as_css");
190  }
191  return ans;
192 }
193 
203 csn *Matrix_as_csn(csn *ans, SEXP x)
204 {
205  static const char *valid[] = {"csn_LU", "csn_QR", ""};
206  int ctype = Matrix_check_class(class_P(x), valid);
207 
208  if (ctype < 0)
209  error(_("invalid class of object to %s"), "Matrix_as_csn");
210  ans->U = Matrix_as_cs(GET_SLOT(x, install("U")));
211  ans->L = Matrix_as_cs(GET_SLOT(x, install("L")));
212  switch(ctype) {
213  case 0:
214  ans->B = (double*) NULL;
215  ans->pinv = INTEGER(GET_SLOT(x, install("Pinv")));
216  break;
217  case 1:
218  ans->B = REAL(GET_SLOT(x, Matrix_betaSym));
219  ans->pinv = (int*) NULL;
220  break;
221  default:
222  error(_("invalid class of object to %s"), "Matrix_as_csn");
223  }
224  return ans;
225 }
226 
239 SEXP Matrix_css_to_SEXP(css *S, char *cl, int dofree, int m, int n)
240 {
241  SEXP ans;
242  static const char *valid[] = {"css_LU", "css_QR", ""};
243  int *nz, ctype = Matrix_check_class(cl, valid);
244 
245  if (ctype < 0)
246  error(_("Inappropriate class cl='%s' in Matrix_css_to_SEXP(S, cl, ..)"),
247  cl);
248  ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cl)));
249  /* allocate and copy common slots */
250  Memcpy(INTEGER(ALLOC_SLOT(ans, install("Q"), INTSXP, n)), S->q, n);
251  nz = INTEGER(ALLOC_SLOT(ans, install("nz"), INTSXP, 3));
252  nz[0] = S->m2; nz[1] = S->lnz; nz[2] = S->unz;
253  switch(ctype) {
254  case 0:
255  break;
256  case 1:
257  Memcpy(INTEGER(ALLOC_SLOT(ans, install("Pinv"), INTSXP, m)),
258  S->pinv, m);
259  Memcpy(INTEGER(ALLOC_SLOT(ans, install("parent"), INTSXP, n)),
260  S->parent, n);
261  Memcpy(INTEGER(ALLOC_SLOT(ans, install("cp"), INTSXP, n)),
262  S->cp, n);
263  break;
264  default:
265  error(_("Inappropriate class cl='%s' in Matrix_css_to_SEXP(S, cl, ..)"),
266  cl);
267  }
268  if (dofree > 0) cs_sfree(S);
269  if (dofree < 0) Free(S);
270  UNPROTECT(1);
271  return ans;
272 }
273 
285 SEXP Matrix_csn_to_SEXP(csn *N, char *cl, int dofree, SEXP dn)
286 {
287  SEXP ans;
288  static const char *valid[] = {"csn_LU", "csn_QR", ""};
289  int ctype = Matrix_check_class(cl, valid), n = (N->U)->n;
290 
291  if (ctype < 0)
292  error(_("Inappropriate class cl='%s' in Matrix_csn_to_SEXP(S, cl, ..)"),
293  cl);
294  ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cl)));
295  /* allocate and copy common slots */
296  /* FIXME: Use the triangular matrix classes for csn_LU */
297  SET_SLOT(ans, install("L"), /* these are free'd later if requested */
298  Matrix_cs_to_SEXP(N->L, "dgCMatrix", 0, dn)); // FIXME: dn
299  SET_SLOT(ans, install("U"),
300  Matrix_cs_to_SEXP(N->U, "dgCMatrix", 0, dn)); // FIXME: dn
301  switch(ctype) {
302  case 0:
303  Memcpy(INTEGER(ALLOC_SLOT(ans, install("Pinv"), INTSXP, n)),
304  N->pinv, n);
305  break;
306  case 1:
307  Memcpy(REAL(ALLOC_SLOT(ans, Matrix_betaSym, REALSXP, n)),
308  N->B, n);
309  break;
310  default:
311  error(_("Inappropriate class cl='%s' in Matrix_csn_to_SEXP(S, cl, ..)"),
312  cl);
313  }
314  if (dofree > 0) cs_nfree(N);
315  if (dofree < 0) {
316  Free(N->L); Free(N->U); Free(N);
317  }
318  UNPROTECT(1);
319  return ans;
320 }
321 
322 #endif /* unused */
SEXP Matrix_DimSym
Definition: Syms.h:2
#define class_P(_x_)
Definition: Mutils.h:178
SEXP Matrix_xSym
Definition: Syms.h:2
csi * pinv
Definition: cs.h:83
csi n
Definition: cs.h:37
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
csi * cp
Definition: cs.h:72
SEXP Matrix_betaSym
Definition: Syms.h:2
double * B
Definition: cs.h:84
static R_INLINE int Matrix_check_class(char *class, const char **valid)
Return the 0-based index of a string match in a vector of strings terminated by an empty string...
Definition: Mutils.h:432
csi m2
Definition: cs.h:74
csi m
Definition: cs.h:36
cholmod_common cl
Definition: chm_common.c:16
SEXP Matrix_uploSym
Definition: Syms.h:2
csi * p
Definition: cs.h:38
cs * U
Definition: cs.h:82
csi * i
Definition: cs.h:39
double lnz
Definition: cs.h:75
Definition: cs.h:79
csi * parent
Definition: cs.h:71
cs * Matrix_as_cs(cs *ans, SEXP x, Rboolean check_Udiag)
Create a cs object with the contents of x.
Definition: cs_utils.c:61
css * cs_sfree(css *S)
Definition: cs.c:1957
double unz
Definition: cs.h:76
static CSP csp_eye(int n)
Create an identity matrix of size n as a cs struct.
Definition: cs_utils.c:31
csi * pinv
Definition: cs.h:69
#define _(String)
Definition: Mutils.h:32
csi nzmax
Definition: cs.h:35
SEXP Matrix_iSym
Definition: Syms.h:2
Definition: cs.h:33
csi * q
Definition: cs.h:70
cs * cs_add(const cs *A, const cs *B, double alpha, double beta)
Definition: cs.c:3
#define diag_P(_x_)
Definition: Mutils.h:175
SEXP Matrix_cs_to_SEXP(cs *a, char *cl, int dofree, SEXP dn)
Copy the contents of a to an appropriate CsparseMatrix object and, optionally, free a or free both a ...
Definition: cs_utils.c:121
csi nz
Definition: cs.h:41
SEXP Matrix_pSym
Definition: Syms.h:2
csn * cs_nfree(csn *N)
Definition: cs.c:1946
cs * cs_spfree(cs *A)
Definition: cs.c:1936
SEXP Matrix_diagSym
Definition: Syms.h:2
double * x
Definition: cs.h:40
static int is_sym(cs *A)
Definition: cs_utils.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
cs * L
Definition: cs.h:81
Definition: cs.h:67
cs * cs_transpose(const cs *A, csi values)
Definition: cs.c:1830
cs * cs_spalloc(csi m, csi n, csi nzmax, csi values, csi triplet)
Definition: cs.c:1907