Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
sparseQR.c
Go to the documentation of this file.
1 #include "sparseQR.h"
2 
3 SEXP sparseQR_validate(SEXP x)
4 {
5  CSP V = AS_CSP__(GET_SLOT(x, Matrix_VSym)),
6  R = AS_CSP__(GET_SLOT(x, install("R")));
7  SEXP beta = GET_SLOT(x, Matrix_betaSym),
8  p = GET_SLOT(x, Matrix_pSym),
9  q = GET_SLOT(x, install("q"));
10  R_CheckStack();
11 
12  if (LENGTH(p) != V->m)
13  return mkString(_("length(p) must match nrow(V)"));
14  if (LENGTH(beta) != V->n)
15  return mkString(_("length(beta) must match ncol(V)"));
16  int lq = LENGTH(q);
17  if (lq && lq != R->n)
18  return mkString(_("length(q) must be zero or ncol(R)"));
19  if (V->n != R->n)
20  return mkString("ncol(V) != ncol(R)");
21  /* FIXME: Check that the permutations are permutations */
22  return ScalarLogical(1);
23 }
24 
36 static
37 void sparseQR_Qmult(cs *V, SEXP dmns, double *beta, int *p, int trans,
38  /* --> */ SEXP ans)
39 {
40  double *y = REAL( GET_SLOT(ans, Matrix_xSym));
41  int *ydims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
42  /* y: contents of a V->m by nrhs, i.e. dim(y) == ydims[0:1], dense matrix
43  * -- Note that V->m = m2 : V may contain "spurious 0 rows" (structural rank deficiency) */
44  int m = V->m, n = V->n;
45  if (ydims[0] != m)
46  error(_("sparseQR_Qmult(): nrow(y) = %d != %d = nrow(V)"), ydims[0], m);
47  double *x; // workspace
48  C_or_Alloca_TO(x, m, double);
49  if (trans) {
50  for (int j = 0; j < ydims[1]; j++) {
51  double *yj = y + j * m;
52  cs_pvec(p, yj, x, m); /* x(0:m-1) = y(p(0:m-1), j) */
53  Memcpy(yj, x, m); /* replace it */
54  for (int k = 0 ; k < n ; k++) /* apply H[1]...H[n] */
55  cs_happly(V, k, beta[k], yj);
56  }
57  } else {
58  for (int j = 0; j < ydims[1]; j++) {
59  double *yj = y + j * m;
60  for (int k = n - 1 ; k >= 0 ; k--) /* apply H[n]...H[1] */
61  cs_happly(V, k, beta[k], yj);
62  cs_ipvec(p, yj, x, m); /* inverse permutation */
63  Memcpy(yj, x, m);
64  }
65  }
66  if(m >= SMALL_4_Alloca) Free(x);
67  if(!isNull(dmns)) { // assign rownames to 'ans' matrix
68  // FIXME? colnames taken from 'y' ?!
69  if(!isNull(VECTOR_ELT(dmns, 0)))
70  SET_VECTOR_ELT(GET_SLOT(ans, Matrix_DimNamesSym), 0,
71  duplicate(VECTOR_ELT(dmns, 0)));
72  }
73 }
74 
75 
84 SEXP sparseQR_qty(SEXP qr, SEXP y, SEXP trans, SEXP keep_dimnames)
85 {
86 
87 //--- will be prepended also to other sparseQR_..() functions below -----------
88 #define INIT_sparseQR_(_DM_NMS_) \
89  SEXP V_ = GET_SLOT(qr, Matrix_VSym); \
90  CSP V = AS_CSP__(V_); \
91  R_CheckStack(); \
92  SEXP ans, aa, dnms = R_NilValue; \
93  if(_DM_NMS_) dnms = GET_SLOT(V_, Matrix_DimNamesSym); \
94  PROTECT_INDEX ipx; \
95  PROTECT_WITH_INDEX(ans = dup_mMatrix_as_dgeMatrix(y), &ipx); \
96  int *ydims = INTEGER(GET_SLOT(ans, Matrix_DimSym)), \
97  m = ydims[0], n = ydims[1], M = V->m, *d_a; \
98  Rboolean rank_def = (m < M); \
99  if(rank_def) { /* must add 0-rows to y, i.e. ans, and remove them *before* return */ \
100  aa = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))); \
101  d_a = INTEGER(GET_SLOT(aa, Matrix_DimSym)); d_a[0] = M; d_a[1] = n; \
102  SEXP dn = GET_SLOT(aa, Matrix_DimNamesSym); \
103  SET_VECTOR_ELT(dn, 1, \
104  duplicate(VECTOR_ELT(GET_SLOT(ans, Matrix_DimNamesSym), 1))); \
105  SET_SLOT(aa, Matrix_DimNamesSym, dn); \
106  double *yy = REAL( GET_SLOT(ans, Matrix_xSym)); /* m * n */ \
107  double *xx = REAL(ALLOC_SLOT(aa, Matrix_xSym, REALSXP, M * n)); \
108  for(int j = 0; j < n; j++) { /* j-th column */ \
109  Memcpy(xx + j*M, yy + j*m, m); /* copy x[ 1:m , j ] := yy[,j] */ \
110  for(int i = m; i < M; i++) xx[i + j*M] = 0.;/* x[(m+1):M, j ] := 0 */ \
111  } \
112  REPROTECT(ans = duplicate(aa), ipx); /* is M x n now */ \
113  }
114 //--- end {INIT_sparseQR_} -----------------------------------------------------
115 
116  INIT_sparseQR_(TRUE)
117  ;
118  sparseQR_Qmult(V, dnms, REAL(GET_SLOT(qr, Matrix_betaSym)),
119  INTEGER(GET_SLOT(qr, Matrix_pSym)), asLogical(trans),
120  ans);
121 #define EXIT_sparseQR_ \
122  /* remove the extra rows from ans */ \
123  d_a[0] = m;/* -> @Dim is ok; @Dimnames (i.e. colnames) still are */ \
124  double *yy = REAL( GET_SLOT(ans, Matrix_xSym)); /* is M x n */ \
125  double *xx = REAL(ALLOC_SLOT(aa, Matrix_xSym, REALSXP, m * n)); \
126  for(int j = 0; j < n; j++) { /* j-th column */ \
127  Memcpy(xx + j*m, yy + j*M, m); /* copy x[ 1:m, j ] := yy[,j] */ \
128  } \
129  ans = duplicate(aa); /* m x n finally */ \
130  UNPROTECT(1) // aa
131 
132 
133  if(rank_def) {
134  warning(_("%s(): structurally rank deficient case: possibly WRONG zeros"),
135  "sparseQR_qty");
137  }
138 
139  UNPROTECT(1);
140  return ans;
141 }
142 
143 // Compute qr.coef(qr, y) := R^{-1} Q' y {modulo row and column permutations}
144 SEXP sparseQR_coef(SEXP qr, SEXP y)
145 {
146  SEXP qslot = GET_SLOT(qr, install("q")), R_ = GET_SLOT(qr, install("R"));
147  CSP R = AS_CSP__(R_);
148 
149  INIT_sparseQR_(FALSE); // (y, ...)
150  // ans := R^{-1} Q' y ==> rownames(ans) := rownames(R^{-1}) = colnames(R) :
151  dnms = duplicate(GET_SLOT(R_, Matrix_DimNamesSym));
152  SET_VECTOR_ELT(dnms, 0, VECTOR_ELT(dnms, 1));
153 
154  /* apply row permutation and multiply by Q' */
155  sparseQR_Qmult(V, dnms, REAL(GET_SLOT(qr, Matrix_betaSym)),
156  INTEGER(GET_SLOT(qr, Matrix_pSym)), /* trans = */ TRUE,
157  ans);
158 
159  double *ax = REAL(GET_SLOT(ans, Matrix_xSym));
160  // FIXME: check n_R, M (= R->m) vs n, m
161  int *q = INTEGER(qslot), lq = LENGTH(qslot), n_R = R->n; // = ncol(R)
162  double *x;
163  if(lq) { C_or_Alloca_TO(x, M, double); }
164  for (int j = 0; j < n; j++) {
165  double *aj = ax + j * M;
166  cs_usolve(R, aj);
167  if (lq) {
168  cs_ipvec(q, aj, x, n_R);
169  Memcpy(aj, x, n_R);
170  }
171  }
172  if(lq && M >= SMALL_4_Alloca) Free(x);
173 
174  if(rank_def) {
175  warning(_("%s(): structurally rank deficient case: possibly WRONG zeros"),
176  "sparseQR_coef");
178  }
179 
180  UNPROTECT(1);
181  return ans;
182 }
183 
186 SEXP sparseQR_resid_fitted(SEXP qr, SEXP y, SEXP want_resid)
187 {
188  int *p = INTEGER(GET_SLOT(qr, Matrix_pSym)),
189  resid = asLogical(want_resid);
190  double *beta = REAL(GET_SLOT(qr, Matrix_betaSym));
191 
192  INIT_sparseQR_(FALSE);
193  // ..... ans should get rownames of 'y' ...
194  /* apply row permutation and multiply by Q' */
195  sparseQR_Qmult(V, dnms, beta, p, /* trans = */ TRUE, ans);
196 
197  double *ax = REAL(GET_SLOT(ans, Matrix_xSym));
198 // FIXME (n,m) := dim(y) vs (N,M) := dim(V) -- ok ??
199  int N = V->n; // M = V->m (in INIT_.. above)
200  for (int j = 0; j < n; j++) {
201  if (resid) // qr.resid(): zero first N rows
202  for (int i = 0; i < N; i++) ax[i + j * M] = 0;
203  else // qr.fitted(): zero last M - N rows
204  for (int i = N; i < M; i++) ax[i + j * M] = 0;
205  }
206  /* multiply by Q and apply inverse row permutation */
207  sparseQR_Qmult(V, dnms, beta, p, /* trans = */ FALSE, ans);
208 
209  if(rank_def) {
210  warning(_("%s(): structurally rank deficient case: possibly WRONG zeros"),
211  "sparseQR_resid_fitted");
213  }
214 
215  UNPROTECT(1);
216  return ans;
217 }
SEXP Matrix_DimSym
Definition: Syms.h:2
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
#define AS_CSP__(x)
Definition: cs_utils.h:13
SEXP Matrix_xSym
Definition: Syms.h:2
csi n
Definition: cs.h:37
SEXP sparseQR_qty(SEXP qr, SEXP y, SEXP trans, SEXP keep_dimnames)
Given a sparse QR decomposition and y, compute Q y or Q'y.
Definition: sparseQR.c:84
SEXP sparseQR_validate(SEXP x)
Definition: sparseQR.c:3
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
SEXP Matrix_betaSym
Definition: Syms.h:2
csi m
Definition: cs.h:36
#define SMALL_4_Alloca
Definition: Mutils.h:47
SEXP sparseQR_resid_fitted(SEXP qr, SEXP y, SEXP want_resid)
Compute qr.resid(qr, y) or qr.fitted(qr, y)
Definition: sparseQR.c:186
SEXP sparseQR_coef(SEXP qr, SEXP y)
Definition: sparseQR.c:144
#define _(String)
Definition: Mutils.h:32
csi cs_usolve(const cs *U, double *x)
Definition: cs.c:1890
Definition: cs.h:33
#define INIT_sparseQR_(_DM_NMS_)
SEXP Matrix_pSym
Definition: Syms.h:2
csi cs_ipvec(const csi *p, const double *b, double *x, csi n)
Definition: cs.c:942
static void sparseQR_Qmult(cs *V, SEXP dmns, double *beta, int *p, int trans, SEXP ans)
Apply Householder transformations and the row permutation P to y.
Definition: sparseQR.c:37
#define EXIT_sparseQR_
SEXP Matrix_VSym
Definition: Syms.h:2
csi cs_happly(const cs *V, csi i, double beta, double *x)
Definition: cs.c:902