Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
dgCMatrix.c
Go to the documentation of this file.
1 #include "dgCMatrix.h"
2 
3 /* for Csparse_transpose() : */
4 #include "Csparse.h"
5 #include "chm_common.h"
6 /* -> Mutils.h / SPQR ... */
7 
8 /* FIXME -- we "forget" about dimnames almost everywhere : */
9 
10 /* for dgCMatrix _and_ lgCMatrix and others (but *not* ngC...) : */
11 SEXP xCMatrix_validate(SEXP x)
12 {
13  /* Almost everything now in Csparse_validate ( ./Csparse.c )
14  * *but* the checking of the 'x' slot : */
15  if (length(GET_SLOT(x, Matrix_iSym)) !=
16  length(GET_SLOT(x, Matrix_xSym)))
17  return mkString(_("lengths of slots 'i' and 'x' must match"));
18 
19  return ScalarLogical(1);
20 }
21 
22 /* for dgRMatrix _and_ lgRMatrix and others (but *not* ngC...) : */
23 SEXP xRMatrix_validate(SEXP x)
24 {
25  /* Almost everything now in Rsparse_validate ( ./Csparse.c )
26  * *but* the checking of the 'x' slot : */
27  if (length(GET_SLOT(x, Matrix_jSym)) !=
28  length(GET_SLOT(x, Matrix_xSym)))
29  return mkString(_("lengths of slots 'j' and 'x' must match"));
30 
31  return ScalarLogical(1);
32 }
33 
34 /* This and the following R_to_CMatrix() lead to memory-not-mapped seg.faults
35  * only with {32bit + R-devel + enable-R-shlib} -- no idea why */
36 SEXP compressed_to_TMatrix(SEXP x, SEXP colP)
37 {
38  int col = asLogical(colP); /* 1 if "C"olumn compressed; 0 if "R"ow */
39  /* however, for Csparse, we now effectively use the cholmod-based
40  * Csparse_to_Tsparse() in ./Csparse.c ; maybe should simply write
41  * an as_cholmod_Rsparse() function and then do "as there" ...*/
42  SEXP indSym = col ? Matrix_iSym : Matrix_jSym,
43  ans, indP = GET_SLOT(x, indSym),
44  pP = GET_SLOT(x, Matrix_pSym);
45  int npt = length(pP) - 1;
46  char *ncl = strdup(class_P(x));
47  static const char *valid[] = { MATRIX_VALID_Csparse, MATRIX_VALID_Rsparse, ""};
48  int ctype = R_check_class_etc(x, valid);
49 
50  if (ctype < 0)
51  error(_("invalid class(x) '%s' in compressed_to_TMatrix(x)"), ncl);
52 
53  /* replace 'C' or 'R' with 'T' :*/
54  ncl[2] = 'T';
55  ans = PROTECT(NEW_OBJECT(MAKE_CLASS(ncl)));
56 
57  slot_dup(ans, x, Matrix_DimSym);
58  if((ctype / 3) % 4 != 2) /* not n..Matrix */
59  slot_dup(ans, x, Matrix_xSym);
60  if(ctype % 3) { /* s(ymmetric) or t(riangular) : */
61  slot_dup(ans, x, Matrix_uploSym);
62  if(ctype % 3 == 2) /* t(riangular) : */
63  slot_dup(ans, x, Matrix_diagSym);
64  }
65  SET_DimNames(ans, x);
66  // possibly asymmetric for symmetricMatrix is ok
67  SET_SLOT(ans, indSym, duplicate(indP));
68  expand_cmprPt(npt, INTEGER(pP),
69  INTEGER(ALLOC_SLOT(ans, col ? Matrix_jSym : Matrix_iSym,
70  INTSXP, length(indP))));
71  free(ncl);
72  UNPROTECT(1);
73  return ans;
74 }
75 
76 SEXP R_to_CMatrix(SEXP x)
77 {
78  SEXP ans, tri = PROTECT(allocVector(LGLSXP, 1));
79  char *ncl = strdup(class_P(x));
80  static const char *valid[] = { MATRIX_VALID_Rsparse, ""};
81  int ctype = R_check_class_etc(x, valid);
82  int *x_dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), *a_dims;
83  PROTECT_INDEX ipx;
84 
85  if (ctype < 0)
86  error(_("invalid class(x) '%s' in R_to_CMatrix(x)"), ncl);
87 
88  /* replace 'R' with 'C' : */
89  ncl[2] = 'C';
90  PROTECT_WITH_INDEX(ans = NEW_OBJECT(MAKE_CLASS(ncl)), &ipx);
91 
92  a_dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
93  /* reversed dim() since we will transpose: */
94  a_dims[0] = x_dims[1];
95  a_dims[1] = x_dims[0];
96 
97  /* triangular: */ LOGICAL(tri)[0] = 0;
98  if((ctype / 3) != 2) /* not n..Matrix */
99  slot_dup(ans, x, Matrix_xSym);
100  if(ctype % 3) { /* s(ymmetric) or t(riangular) : */
101  SET_SLOT(ans, Matrix_uploSym,
102  mkString((*uplo_P(x) == 'U') ? "L" : "U"));
103  if(ctype % 3 == 2) { /* t(riangular) : */
104  LOGICAL(tri)[0] = 1;
105  slot_dup(ans, x, Matrix_diagSym);
106  }
107  }
108  SET_SLOT(ans, Matrix_iSym, duplicate(GET_SLOT(x, Matrix_jSym)));
109  slot_dup(ans, x, Matrix_pSym);
110  REPROTECT(ans = Csparse_transpose(ans, tri), ipx);
111  SET_DimNames(ans, x);
112  // possibly asymmetric for symmetricMatrix is ok
113  free(ncl);
114  UNPROTECT(2);
115  return ans;
116 }
117 
122 SEXP compressed_non_0_ij(SEXP x, SEXP colP)
123 {
124  int col = asLogical(colP); /* 1 if "C"olumn compressed; 0 if "R"ow */
125  SEXP ans, indSym = col ? Matrix_iSym : Matrix_jSym;
126  SEXP indP = GET_SLOT(x, indSym),
127  pP = GET_SLOT(x, Matrix_pSym);
128  int i, *ij;
129  int nouter = INTEGER(GET_SLOT(x, Matrix_DimSym))[col ? 1 : 0],
130  n_el = INTEGER(pP)[nouter]; /* is only == length(indP), if the
131  inner slot is not over-allocated */
132 
133  ij = INTEGER(ans = PROTECT(allocMatrix(INTSXP, n_el, 2)));
134  /* expand the compressed margin to 'i' or 'j' : */
135  expand_cmprPt(nouter, INTEGER(pP), &ij[col ? n_el : 0]);
136  /* and copy the other one: */
137  if (col)
138  for(i = 0; i < n_el; i++)
139  ij[i] = INTEGER(indP)[i];
140  else /* row compressed */
141  for(i = 0; i < n_el; i++)
142  ij[i + n_el] = INTEGER(indP)[i];
143 
144  UNPROTECT(1);
145  return ans;
146 }
147 
148 #if 0 /* unused */
149 SEXP dgCMatrix_lusol(SEXP x, SEXP y)
150 {
151  SEXP ycp = PROTECT((TYPEOF(y) == REALSXP) ?
152  duplicate(y) : coerceVector(y, REALSXP));
153  CSP xc = AS_CSP__(x);
154  R_CheckStack();
155 
156  if (xc->m != xc->n || xc->m <= 0)
157  error(_("dgCMatrix_lusol requires a square, non-empty matrix"));
158  if (LENGTH(ycp) != xc->m)
159  error(_("Dimensions of system to be solved are inconsistent"));
160  if (!cs_lusol(/*order*/ 1, xc, REAL(ycp), /*tol*/ 1e-7))
161  error(_("cs_lusol failed"));
162 
163  UNPROTECT(1);
164  return ycp;
165 }
166 #endif
167 
168 // called from package MatrixModels's R code
169 SEXP dgCMatrix_qrsol(SEXP x, SEXP y, SEXP ord)
170 {
171  /* FIXME: extend this to work in multivariate case, i.e. y a matrix with > 1 column ! */
172  SEXP ycp = PROTECT((TYPEOF(y) == REALSXP) ?
173  duplicate(y) : coerceVector(y, REALSXP));
174  CSP xc = AS_CSP(x); /* <--> x may be dgC* or dtC* */
175  int order = asInteger(ord);
176 #ifdef _not_yet_do_FIXME__
177  const char *nms[] = {"L", "coef", "Xty", "resid", ""};
178  SEXP ans = PROTECT(Rf_mkNamed(VECSXP, nms));
179 #endif
180  R_CheckStack();
181 
182  if (order < 0 || order > 3)
183  error(_("dgCMatrix_qrsol(., order) needs order in {0,..,3}"));
184  /* --> cs_amd() --- order 0: natural, 1: Chol, 2: LU, 3: QR */
185  if (LENGTH(ycp) != xc->m)
186  error(_("Dimensions of system to be solved are inconsistent"));
187  /* FIXME? Note that qr_sol() would allow *under-determined systems;
188  * In general, we'd need LENGTH(ycp) = max(n,m)
189  * FIXME also: multivariate y (see above)
190  */
191  if (xc->m < xc->n || xc->n <= 0)
192  error(_("dgCMatrix_qrsol(<%d x %d>-matrix) requires a 'tall' rectangular matrix"),
193  xc->m, xc->n);
194 
195  /* cs_qrsol(): Tim Davis (2006) .. "8.2 Using a QR factorization", p.136f , calling
196  * ------- cs_sqr(order, ..), see p.76 */
197  /* MM: FIXME: write our *OWN* version of - the first case (m >= n) - of cs_qrsol()
198  * --------- which will (1) work with a *multivariate* y
199  * (2) compute coefficients properly, not overwriting RHS
200  */
201  if (!cs_qrsol(order, xc, REAL(ycp)))
202  /* return value really is 0 or 1 - no more info there */
203  error(_("cs_qrsol() failed inside dgCMatrix_qrsol()"));
204 
205  /* Solution is only in the first part of ycp -- cut its length back to n : */
206  ycp = lengthgets(ycp, (R_xlen_t) xc->n);
207 
208  UNPROTECT(1);
209  return ycp;
210 }
211 
212 // Modified version of Tim Davis's cs_qr_mex.c file for MATLAB (in CSparse)
213 // Usage: [V,beta,p,R,q] = cs_qr(A) ;
214 SEXP dgCMatrix_QR(SEXP Ap, SEXP order, SEXP keep_dimnames)
215 {
216  CSP A = AS_CSP__(Ap), D;
217  int io = INTEGER(order)[0];
218  Rboolean verbose = (io < 0);// verbose=TRUE, encoded with negative 'order'
219  int m0 = A->m, m = m0, n = A->n, ord = asLogical(order) ? 3 : 0, *p;
220  R_CheckStack();
221 
222  if (m < n) error(_("A must have #{rows} >= #{columns}")) ;
223  SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sparseQR")));
224  int *dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
225  dims[0] = m; dims[1] = n;
226  css *S = cs_sqr(ord, A, 1); /* symbolic QR ordering & analysis*/
227  if (!S) error(_("cs_sqr failed"));
228  int keep_dimnms = asLogical(keep_dimnames);
229  if(keep_dimnms == NA_LOGICAL) { keep_dimnms = TRUE;
230  warning(_("dgcMatrix_QR(*, keep_dimnames = NA): NA taken as TRUE"));
231  }
232  if(verbose && S->m2 > m) // in ./cs.h , m2 := # of rows for QR, after adding fictitious rows
233  Rprintf("Symbolic QR(): Matrix structurally rank deficient (m2-m = %d)\n",
234  S->m2 - m);
235  csn *N = cs_qr(A, S); /* numeric QR factorization */
236  if (!N) error(_("cs_qr failed")) ;
237  cs_dropzeros(N->L); /* drop zeros from V and sort */
238  D = cs_transpose(N->L, 1); cs_spfree(N->L);
239  N->L = cs_transpose(D, 1); cs_spfree(D);
240  cs_dropzeros(N->U); /* drop zeros from R and sort */
241  D = cs_transpose(N->U, 1); cs_spfree(N->U) ;
242  N->U = cs_transpose(D, 1); cs_spfree(D);
243  m = N->L->m; /* m may be larger now */
244  // MM: m := S->m2 also counting the ficticious rows (Tim Davis, p.72, 74f)
245  p = cs_pinv(S->pinv, m); /* p = pinv' */
246  SEXP dn = R_NilValue; Rboolean do_dn = FALSE;
247  if(keep_dimnms) {
248  dn = GET_SLOT(Ap, Matrix_DimNamesSym);
249  do_dn = !isNull(VECTOR_ELT(dn, 0)) && m == m0;
250  // FIXME? also deal with case m > m0 ?
251  if(do_dn) { // keep rownames
252  dn = PROTECT(duplicate(dn));
253  SET_VECTOR_ELT(dn, 1, R_NilValue);
254  } else dn = R_NilValue;
255  }
256  SET_SLOT(ans, Matrix_VSym, Matrix_cs_to_SEXP(N->L, "dgCMatrix", 0, dn)); // "V"
257  Memcpy(REAL(ALLOC_SLOT(ans, Matrix_betaSym, REALSXP, n)), N->B, n);
258  Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, m)), p, m);
259  if(do_dn) {
260  UNPROTECT(1); // dn
261  dn = R_NilValue; do_dn = FALSE;
262  }
263  if (ord) {
264  Memcpy(INTEGER(ALLOC_SLOT(ans, install("q"), INTSXP, n)), S->q, n);
265  if(keep_dimnms) {
266  dn = GET_SLOT(Ap, Matrix_DimNamesSym);
267  do_dn = !isNull(VECTOR_ELT(dn, 1));
268  if(do_dn) {
269  dn = PROTECT(duplicate(dn));
270  // permute colnames by S->q : cn <- cn[ S->q ] :
271  SEXP cns = PROTECT(duplicate(VECTOR_ELT(dn, 1)));
272  for(int j=0; j < n; j++)
273  SET_STRING_ELT(VECTOR_ELT(dn, 1), j, STRING_ELT(cns, S->q[j]));
274  UNPROTECT(1);
275  SET_VECTOR_ELT(dn, 0, R_NilValue);
276  } else dn = R_NilValue;
277  }
278  } else
279  ALLOC_SLOT(ans, install("q"), INTSXP, 0);
280  SET_SLOT(ans, install("R"), Matrix_cs_to_SEXP(N->U, "dgCMatrix", 0, dn));
281  if(do_dn) UNPROTECT(1); // dn
282  cs_nfree(N);
283  cs_sfree(S);
284  cs_free(p);
285  UNPROTECT(1);
286  return ans;
287 }
288 
289 #ifdef Matrix_with_SPQR
290 
310 SEXP dgCMatrix_SPQR(SEXP Ap, SEXP ordering, SEXP econ, SEXP tol)
311 {
312 /* SEXP ans = PROTECT(allocVector(VECSXP, 4)); */
313  SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("SPQR")));
314 
315  CHM_SP A = AS_CHM_SP(Ap), Q, R;
316  SuiteSparse_long *E, rank;/* not always = int FIXME (Windows_64 ?) */
317 
318  if ((rank = SuiteSparseQR_C_QR(asInteger(ordering),
319  asReal(tol),/* originally had SPQR_DEFAULT_TOL */
320  (SuiteSparse_long)asInteger(econ),/* originally had 0 */
321  A, &Q, &R, &E, &cl)) == -1)
322  error(_("SuiteSparseQR_C_QR returned an error code"));
323 
324  slot_dup(ans, Ap, Matrix_DimSym);
325 /* SET_VECTOR_ELT(ans, 0, */
326 /* chm_sparse_to_SEXP(Q, 0, 0, 0, "", R_NilValue)); */
327  SET_SLOT(ans, install("Q"),
328  chm_sparse_to_SEXP(Q, 0, 0, 0, "", R_NilValue));
329 
330  /* Also gives a dgCMatrix (not a dtC* *triangular*) :
331  * may make sense if to be used in the "spqr_solve" routines .. ?? */
332 /* SET_VECTOR_ELT(ans, 1, */
333 /* chm_sparse_to_SEXP(R, 0, 0, 0, "", R_NilValue)); */
334  SET_SLOT(ans, install("R"),
335  chm_sparse_to_SEXP(R, 0, 0, 0, "", R_NilValue));
336  cholmod_free_sparse(&Al, &cl);
337  cholmod_free_sparse(&R, &cl);
338  cholmod_free_sparse(&Q, &cl);
339  if (E) {
340  int *Er;
341  SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, A->ncol));
342  Er = INTEGER(VECTOR_ELT(ans, 2));
343  for (int i = 0; i < A->ncol; i++) Er[i] = (int) E[i];
344  Free(E);
345  } else SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, 0));
346  SET_VECTOR_ELT(ans, 3, ScalarInteger((int)rank));
347  UNPROTECT(1);
348  return ans;
349 }
350 #endif
351 /* Matrix_with_SPQR */
352 
353 /* Modified version of Tim Davis's cs_lu_mex.c file for MATLAB */
354 void install_lu(SEXP Ap, int order, double tol, Rboolean err_sing, Rboolean keep_dimnms)
355 {
356  // (order, tol) == (1, 1) by default, when called from R.
357  SEXP ans;
358  css *S;
359  csn *N;
360  int n, *p, *dims;
361  CSP A = AS_CSP__(Ap), D;
362  R_CheckStack();
363 
364  n = A->n;
365  if (A->m != n)
366  error(_("LU decomposition applies only to square matrices"));
367  if (order) { /* not using natural order */
368  order = (tol == 1) ? 2 /* amd(S'*S) w/dense rows or I */
369  : 1; /* amd (A+A'), or natural */
370  }
371  S = cs_sqr(order, A, /*qr = */ 0); /* symbolic ordering */
372  N = cs_lu(A, S, tol); /* numeric factorization */
373  if (!N) {
374  if(err_sing)
375  error(_("cs_lu(A) failed: near-singular A (or out of memory)"));
376  else {
377  /* No warning: The useR should be careful :
378  * Put NA into "LU" factor */
379  set_factors(Ap, ScalarLogical(NA_LOGICAL), "LU");
380  return;
381  }
382  }
383  cs_dropzeros(N->L); /* drop zeros from L and sort it */
384  D = cs_transpose(N->L, 1);
385  cs_spfree(N->L);
386  N->L = cs_transpose(D, 1);
387  cs_spfree(D);
388  cs_dropzeros(N->U); /* drop zeros from U and sort it */
389  D = cs_transpose(N->U, 1);
390  cs_spfree(N->U);
391  N->U = cs_transpose(D, 1);
392  cs_spfree(D);
393  p = cs_pinv(N->pinv, n); /* p=pinv' */
394  ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sparseLU")));
395  dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
396  dims[0] = n; dims[1] = n;
397  SEXP dn; Rboolean do_dn = FALSE;
398  if(keep_dimnms) {
399  dn = GET_SLOT(Ap, Matrix_DimNamesSym);
400  do_dn = !isNull(VECTOR_ELT(dn, 0));
401  if(do_dn) {
402  dn = PROTECT(duplicate(dn));
403  // permute rownames by p : rn <- rn[ p ] :
404  SEXP rn = PROTECT(duplicate(VECTOR_ELT(dn, 0)));
405  for(int i=0; i < n; i++)
406  SET_STRING_ELT(VECTOR_ELT(dn, 0), i, STRING_ELT(rn, p[i]));
407  UNPROTECT(1); // rn
408  SET_VECTOR_ELT(dn, 1, R_NilValue); // colnames(.) := NULL
409  }
410  }
411  SET_SLOT(ans, install("L"),
412  Matrix_cs_to_SEXP(N->L, "dtCMatrix", 0, do_dn ? dn : R_NilValue));
413 
414  if(keep_dimnms) {
415  if(do_dn) {
416  UNPROTECT(1); // dn
417  dn = GET_SLOT(Ap, Matrix_DimNamesSym);
418  }
419  do_dn = !isNull(VECTOR_ELT(dn, 1));
420  if(do_dn) {
421  dn = PROTECT(duplicate(dn));
422  if(order) { // permute colnames by S->q : cn <- cn[ S->q ] :
423  SEXP cn = PROTECT(duplicate(VECTOR_ELT(dn, 1)));
424  for(int j=0; j < n; j++)
425  SET_STRING_ELT(VECTOR_ELT(dn, 1), j, STRING_ELT(cn, S->q[j]));
426  UNPROTECT(1); // cn
427  }
428  SET_VECTOR_ELT(dn, 0, R_NilValue); // rownames(.) := NULL
429  }
430  }
431  SET_SLOT(ans, install("U"),
432  Matrix_cs_to_SEXP(N->U, "dtCMatrix", 0, do_dn ? dn : R_NilValue));
433  if(do_dn) UNPROTECT(1); // dn
434  Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, /* "p" */
435  INTSXP, n)), p, n);
436  if (order)
437  Memcpy(INTEGER(ALLOC_SLOT(ans, install("q"),
438  INTSXP, n)), S->q, n);
439  cs_nfree(N);
440  cs_sfree(S);
441  cs_free(p);
442  UNPROTECT(1);
443  set_factors(Ap, ans, "LU");
444 }
445 
446 SEXP dgCMatrix_LU(SEXP Ap, SEXP orderp, SEXP tolp, SEXP error_on_sing, SEXP keep_dimnames)
447 {
448  SEXP ans;
449  Rboolean err_sing = asLogical(error_on_sing);
450  /* FIXME: dgCMatrix_LU should check ans for consistency in
451  * permutation type with the requested value - Should have two
452  * classes or two different names in the factors list for LU with
453  * permuted columns or not.
454  * OTOH, currently (order, tol) === (1, 1) always.
455  * It is true that length(LU@q) does flag the order argument.
456  */
457  if (!isNull(ans = get_factors(Ap, "LU")))
458  return ans;
459  int keep_dimnms = asLogical(keep_dimnames);
460  if(keep_dimnms == NA_LOGICAL) { keep_dimnms = TRUE;
461  warning(_("dgcMatrix_LU(*, keep_dimnames = NA): NA taken as TRUE"));
462  }
463  install_lu(Ap, asInteger(orderp), asReal(tolp), err_sing, keep_dimnms);
464  return get_factors(Ap, "LU");
465 }
466 
467 SEXP dgCMatrix_matrix_solve(SEXP Ap, SEXP b, SEXP give_sparse)
468 // FIXME: add 'keep_dimnames' as argument
469 {
470  Rboolean sparse = asLogical(give_sparse);
471  if(sparse) {
472  // FIXME: implement this
473  error(_("dgCMatrix_matrix_solve(.., sparse=TRUE) not yet implemented"));
474 
475  /* Idea: in the for(j = 0; j < nrhs ..) loop below, build the *sparse* result matrix
476  * ----- *column* wise -- which is perfect for dgCMatrix
477  * --> build (i,p,x) slots "increasingly" [well, allocate in batches ..]
478  *
479  * --> maybe first a protoype in R
480  */
481 
482  }
483  SEXP ans = PROTECT(dup_mMatrix_as_dgeMatrix(b)),
484  lu, qslot;
485  CSP L, U;
486  int *bdims = INTEGER(GET_SLOT(ans, Matrix_DimSym)), *p, *q;
487  int j, n = bdims[0], nrhs = bdims[1];
488  double *x, *ax = REAL(GET_SLOT(ans, Matrix_xSym));
489  C_or_Alloca_TO(x, n, double);
490 
491  if (isNull(lu = get_factors(Ap, "LU"))) {
492  install_lu(Ap, /* order = */ 1, /* tol = */ 1.0, /* err_sing = */ TRUE,
493  /* keep_dimnames = */ TRUE);
494  lu = get_factors(Ap, "LU");
495  }
496  qslot = GET_SLOT(lu, install("q"));
497  L = AS_CSP__(GET_SLOT(lu, install("L")));
498  U = AS_CSP__(GET_SLOT(lu, install("U")));
499  R_CheckStack();
500  if (U->n != n)
501  error(_("Dimensions of system to be solved are inconsistent"));
502  if(nrhs >= 1 && n >= 1) {
503  p = INTEGER(GET_SLOT(lu, Matrix_pSym));
504  q = LENGTH(qslot) ? INTEGER(qslot) : (int *) NULL;
505 
506  for (j = 0; j < nrhs; j++) {
507  cs_pvec(p, ax + j * n, x, n); /* x = b(p) */
508  cs_lsolve(L, x); /* x = L\x */
509  cs_usolve(U, x); /* x = U\x */
510  if (q) /* r(q) = x , hence
511  r = Q' U{^-1} L{^-1} P b = A^{-1} b */
512  cs_ipvec(q, x, ax + j * n, n);
513  else
514  Memcpy(ax + j * n, x, n);
515  }
516  }
517  if(n >= SMALL_4_Alloca) Free(x);
518  UNPROTECT(1);
519  return ans;
520 }
521 
522 // called from package MatrixModels's R code:
523 SEXP dgCMatrix_cholsol(SEXP x, SEXP y)
524 {
525  /* Solve Sparse Least Squares X %*% beta ~= y with dense RHS y,
526  * where X = t(x) i.e. we pass x = t(X) as argument,
527  * via "Cholesky(X'X)" .. well not really:
528  * cholmod_factorize("x", ..) finds L in X'X = L'L directly */
529  CHM_SP cx = AS_CHM_SP(x);
530  /* FIXME: extend this to work in multivariate case, i.e. y a matrix with > 1 column ! */
531  CHM_DN cy = AS_CHM_DN(coerceVector(y, REALSXP)), rhs, cAns, resid;
532  CHM_FR L;
533  int n = cx->ncol;/* #{obs.} {x = t(X) !} */
534  double one[] = {1,0}, zero[] = {0,0}, neg1[] = {-1,0};
535  const char *nms[] = {"L", "coef", "Xty", "resid", ""};
536  SEXP ans = PROTECT(Rf_mkNamed(VECSXP, nms));
537  R_CheckStack();
538 
539  if (n < cx->nrow || n <= 0)
540  error(_("dgCMatrix_cholsol requires a 'short, wide' rectangular matrix"));
541  if (cy->nrow != n)
542  error(_("Dimensions of system to be solved are inconsistent"));
543  rhs = cholmod_allocate_dense(cx->nrow, 1, cx->nrow, CHOLMOD_REAL, &c);
544  /* cholmod_sdmult(A, transp, alpha, beta, X, Y, &c):
545  * Y := alpha*(A*X) + beta*Y or alpha*(A'*X) + beta*Y ;
546  * here: rhs := 1 * x %*% y + 0 = x %*% y = X'y */
547  if (!(cholmod_sdmult(cx, 0 /* trans */, one, zero, cy, rhs, &c)))
548  error(_("cholmod_sdmult error (rhs)"));
549  L = cholmod_analyze(cx, &c);
550  if (!cholmod_factorize(cx, L, &c))
551  error(_("cholmod_factorize failed: status %d, minor %d from ncol %d"),
552  c.status, L->minor, L->n);
553 /* FIXME: Do this in stages so an "effects" vector can be calculated */
554  if (!(cAns = cholmod_solve(CHOLMOD_A, L, rhs, &c)))
555  error(_("cholmod_solve (CHOLMOD_A) failed: status %d, minor %d from ncol %d"),
556  c.status, L->minor, L->n);
557  /* L : */
558  SET_VECTOR_ELT(ans, 0, chm_factor_to_SEXP(L, 0));
559  /* coef : */
560  SET_VECTOR_ELT(ans, 1, allocVector(REALSXP, cx->nrow));
561  Memcpy(REAL(VECTOR_ELT(ans, 1)), (double*)(cAns->x), cx->nrow);
562  /* X'y : */
563 /* FIXME: Change this when the "effects" vector is available */
564  SET_VECTOR_ELT(ans, 2, allocVector(REALSXP, cx->nrow));
565  Memcpy(REAL(VECTOR_ELT(ans, 2)), (double*)(rhs->x), cx->nrow);
566  /* resid := y */
567  resid = cholmod_copy_dense(cy, &c);
568  /* cholmod_sdmult(A, transp, alp, bet, X, Y, &c):
569  * Y := alp*(A*X) + bet*Y or alp*(A'*X) + beta*Y ;
570  * here: resid := -1 * x' %*% coef + 1 * y = y - X %*% coef */
571  if (!(cholmod_sdmult(cx, 1/* trans */, neg1, one, cAns, resid, &c)))
572  error(_("cholmod_sdmult error (resid)"));
573  /* FIXME: for multivariate case, i.e. resid *matrix* with > 1 column ! */
574  SET_VECTOR_ELT(ans, 3, allocVector(REALSXP, n));
575  Memcpy(REAL(VECTOR_ELT(ans, 3)), (double*)(resid->x), n);
576 
577  cholmod_free_factor(&L, &c);
578  cholmod_free_dense(&rhs, &c);
579  cholmod_free_dense(&cAns, &c);
580  UNPROTECT(1);
581  return ans;
582 }
583 
584 
585 /* Define all of
586  * dgCMatrix_colSums(....)
587  * igCMatrix_colSums(....)
588  * lgCMatrix_colSums_d(....)
589  * lgCMatrix_colSums_i(....)
590  * ngCMatrix_colSums_d(....)
591  * ngCMatrix_colSums_i(....)
592  */
593 #define _dgC_
594 #include "t_gCMatrix_colSums.c"
595 
596 #define _igC_
597 #include "t_gCMatrix_colSums.c"
598 
599 #define _lgC_
600 #include "t_gCMatrix_colSums.c"
601 
602 #define _ngC_
603 #include "t_gCMatrix_colSums.c"
604 
605 #define _lgC_mn
606 #include "t_gCMatrix_colSums.c"
607 
608 #define _ngC_mn
609 #include "t_gCMatrix_colSums.c"
610 
611 
612 SEXP lgCMatrix_colSums(SEXP x, SEXP NArm, SEXP spRes, SEXP trans, SEXP means)
613 {
614  if(asLogical(means)) /* ==> result will be "double" / "dsparseVector" */
615  return lgCMatrix_colSums_d(x, NArm, spRes, trans, means);
616  else
617  return lgCMatrix_colSums_i(x, NArm, spRes, trans, means);
618 }
619 
620 SEXP ngCMatrix_colSums(SEXP x, SEXP NArm, SEXP spRes, SEXP trans, SEXP means)
621 {
622  if(asLogical(means)) /* ==> result will be "double" / "dsparseVector" */
623  return ngCMatrix_colSums_d(x, NArm, spRes, trans, means);
624  else
625  return ngCMatrix_colSums_i(x, NArm, spRes, trans, means);
626 }
SEXP Matrix_DimSym
Definition: Syms.h:2
csn * cs_lu(const cs *A, const css *S, double tol)
Definition: cs.c:1019
csi cs_qrsol(csi order, const cs *A, double *b)
Definition: cs.c:1477
#define class_P(_x_)
Definition: Mutils.h:178
csi cs_pvec(const csi *p, const double *b, double *x, csi n)
Definition: cs.c:1397
#define C_or_Alloca_TO(_VAR_, _N_, _TYPE_)
Definition: Mutils.h:50
cholmod_sparse * CHM_SP
Definition: chm_common.h:25
csn * cs_qr(const cs *A, const css *S)
Definition: cs.c:1405
#define AS_CSP__(x)
Definition: cs_utils.h:13
SEXP Matrix_xSym
Definition: Syms.h:2
#define slot_dup(dest, src, sym)
Definition: Mutils.h:149
SEXP R_to_CMatrix(SEXP x)
Definition: dgCMatrix.c:76
csi * pinv
Definition: cs.h:83
csi n
Definition: cs.h:37
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
SEXP Matrix_betaSym
Definition: Syms.h:2
double * B
Definition: cs.h:84
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
csi m2
Definition: cs.h:74
#define MATRIX_VALID_Rsparse
Definition: Mutils.h:396
csi m
Definition: cs.h:36
cholmod_common cl
Definition: chm_common.c:16
SEXP Matrix_uploSym
Definition: Syms.h:2
SEXP dup_mMatrix_as_dgeMatrix(SEXP A)
Definition: Mutils.c:852
#define AS_CSP(x)
Definition: cs_utils.h:12
void install_lu(SEXP Ap, int order, double tol, Rboolean err_sing, Rboolean keep_dimnms)
Definition: dgCMatrix.c:354
void * cs_free(void *p)
Definition: cs.c:1148
cs * U
Definition: cs.h:82
SEXP compressed_non_0_ij(SEXP x, SEXP colP)
Return a 2 column matrix '' cbind(i, j) '' of 0-origin index vectors (i,j) which entirely correspond ...
Definition: dgCMatrix.c:122
csi cs_lusol(csi order, const cs *A, double *b, double tol)
Definition: cs.c:1104
SEXP Matrix_jSym
Definition: Syms.h:2
#define SMALL_4_Alloca
Definition: Mutils.h:47
Definition: cs.h:79
csi * cs_pinv(csi const *p, csi n)
Definition: cs.c:1326
SEXP set_factors(SEXP obj, SEXP val, char *nm)
Caches 'val' in the 'factors' slot of obj, i.e.
Definition: Mutils.c:130
SEXP get_factors(SEXP obj, char *nm)
Definition: Mutils.c:106
css * cs_sfree(css *S)
Definition: cs.c:1957
SEXP xCMatrix_validate(SEXP x)
Definition: dgCMatrix.c:11
SEXP xRMatrix_validate(SEXP x)
Definition: dgCMatrix.c:23
cholmod_factor * CHM_FR
Definition: chm_common.h:23
SEXP lgCMatrix_colSums(SEXP x, SEXP NArm, SEXP spRes, SEXP trans, SEXP means)
Definition: dgCMatrix.c:612
static R_INLINE int * expand_cmprPt(int ncol, const int mp[], int mj[])
Expand compressed pointers in the array mp into a full set of indices in the array mj...
Definition: Mutils.h:259
csi cs_dropzeros(cs *A)
Definition: cs.c:761
csi * pinv
Definition: cs.h:69
#define uplo_P(_x_)
Definition: Mutils.h:174
#define _(String)
Definition: Mutils.h:32
csi cs_usolve(const cs *U, double *x)
Definition: cs.c:1890
cholmod_dense * CHM_DN
Definition: chm_common.h:21
SEXP Matrix_iSym
Definition: Syms.h:2
SEXP dgCMatrix_qrsol(SEXP x, SEXP y, SEXP ord)
Definition: dgCMatrix.c:169
Definition: cs.h:33
csi * q
Definition: cs.h:70
SEXP chm_factor_to_SEXP(CHM_FR f, int dofree)
Copy the contents of f to an appropriate CHMfactor object and, optionally, free f or free both f and ...
Definition: chm_common.c:1106
SEXP ngCMatrix_colSums(SEXP x, SEXP NArm, SEXP spRes, SEXP trans, SEXP means)
Definition: dgCMatrix.c:620
#define AS_CHM_DN(x)
Definition: chm_common.h:43
SEXP dgCMatrix_cholsol(SEXP x, SEXP y)
Definition: dgCMatrix.c:523
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
SEXP Matrix_pSym
Definition: Syms.h:2
SEXP dgCMatrix_QR(SEXP Ap, SEXP order, SEXP keep_dimnames)
Definition: dgCMatrix.c:214
csi cs_ipvec(const csi *p, const double *b, double *x, csi n)
Definition: cs.c:942
SEXP compressed_to_TMatrix(SEXP x, SEXP colP)
Definition: dgCMatrix.c:36
csn * cs_nfree(csn *N)
Definition: cs.c:1946
cs * cs_spfree(cs *A)
Definition: cs.c:1936
SEXP Csparse_transpose(SEXP x, SEXP tri)
Definition: Csparse.c:402
SEXP Matrix_diagSym
Definition: Syms.h:2
SEXP dgCMatrix_LU(SEXP Ap, SEXP orderp, SEXP tolp, SEXP error_on_sing, SEXP keep_dimnames)
Definition: dgCMatrix.c:446
cholmod_common c
Definition: chm_common.c:15
csi cs_lsolve(const cs *L, double *x)
Definition: cs.c:985
static R_INLINE void SET_DimNames(SEXP dest, SEXP src)
Definition: Mutils.h:161
#define MATRIX_VALID_Csparse
Definition: Mutils.h:384
css * cs_sqr(csi order, const cs *A, csi qr)
Definition: cs.c:1740
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
cs * L
Definition: cs.h:81
SEXP Matrix_VSym
Definition: Syms.h:2
Definition: cs.h:67
cs * cs_transpose(const cs *A, csi values)
Definition: cs.c:1830
SEXP dgCMatrix_matrix_solve(SEXP Ap, SEXP b, SEXP give_sparse)
Definition: dgCMatrix.c:467