Matrix r5059
Loading...
Searching...
No Matches
Csparse.c
Go to the documentation of this file.
1#include "Mdefines.h"
2#include "cs-etc.h"
3#include "cholmod-etc.h"
4
5/* NB: mostly parallel to CsparseMatrix_validate in ./validity.c */
6SEXP checkpi(SEXP dim, SEXP p, SEXP i)
7{
8 int m = INTEGER(dim)[0], n = INTEGER(dim)[1];
9 if (TYPEOF(p) != INTSXP)
10 return errorString(_("'%s' slot is not of type \"%s\""),
11 "p", "integer");
12 if (XLENGTH(p) - 1 != n)
13 return errorString(_("'%s' slot does not have length %s"),
14 "p", "Dim[2]+1");
15 int *pp = INTEGER(p);
16 if (pp[0] != 0)
17 return errorString(_("first element of '%s' slot is not 0"),
18 "p");
19 int j;
20 for (j = 1; j <= n; ++j) {
21 if (pp[j] == NA_INTEGER)
22 return errorString(_("'%s' slot contains NA"),
23 "p");
24 if (pp[j] < pp[j - 1])
25 return errorString(_("'%s' slot is not nondecreasing"),
26 "p");
27 if (pp[j] - pp[j - 1] > m)
28 return errorString(_("first differences of '%s' slot exceed %s"),
29 "p", "Dim[1]");
30 }
31 if (TYPEOF(i) != INTSXP)
32 return errorString(_("'%s' slot is not of type \"%s\""),
33 "i", "integer");
34 if (XLENGTH(i) < pp[n])
35 return errorString(_("'%s' slot has length less than %s"),
36 "i", "p[length(p)]");
37 int *pi = INTEGER(i), k, kend, ik, i0, sorted = 1;
38 for (j = 1, k = 0; j <= n; ++j) {
39 kend = pp[j];
40 i0 = -1;
41 while (k < kend) {
42 ik = pi[k];
43 if (ik == NA_INTEGER)
44 return errorString(_("'%s' slot contains NA"),
45 "i");
46 if (ik < 0 || ik >= m)
47 return errorString(_("'%s' slot has elements not in {%s}"),
48 "i", "0,...,Dim[1]-1");
49 if (ik < i0)
50 sorted = 0;
51 else if (ik == i0)
52 return errorString(_("'%s' slot is not increasing within columns after sorting"),
53 "i");
54 i0 = ik;
55 ++k;
56 }
57 }
58 SEXP ans = Rf_allocVector(LGLSXP, 1);
59 LOGICAL(ans)[0] = sorted;
60 return ans;
61}
62
63/* .validateCsparse(x, sort.if.needed = TRUE) */
64SEXP CsparseMatrix_validate_maybe_sorting(SEXP x)
65{
66 SEXP dim = PROTECT(GET_SLOT(x, Matrix_DimSym)),
67 p = PROTECT(GET_SLOT(x, Matrix_pSym)),
68 i = PROTECT(GET_SLOT(x, Matrix_iSym)),
69 cpi = checkpi(dim, p, i);
70 if (TYPEOF(cpi) == LGLSXP && !LOGICAL(cpi)[0]) {
71 cholmod_sparse *A = M2CHS(x, 1);
72 A->sorted = 0;
73 if (!cholmod_sort(A, &c))
74 Rf_error(_("'%s' failed"), "cholmod_sort");
75 cpi = checkpi(dim, p, i);
76 }
77 UNPROTECT(3); /* i, p, dim */
78 return cpi;
79}
80
81/* TODO: support NCOL(b) > 1 */
82SEXP dgCMatrix_lusol(SEXP a, SEXP b)
83{
84 Matrix_cs *A = M2CXS(a, 1);
85 CXSPARSE_XTYPE_SET(CXSPARSE_REAL);
86 PROTECT(b = (TYPEOF(b) == REALSXP) ?
87 Rf_duplicate(b) : Rf_coerceVector(b, REALSXP));
88 if (A->m != A->n || A->m <= 0)
89 Rf_error(_("'%s' is empty or not square"), "a");
90 if (LENGTH(b) != A->m)
91 Rf_error(_("dimensions of '%s' and '%s' are inconsistent"), "a", "b");
92 if (!Matrix_cs_lusol(1, A, REAL(b), 1e-07))
93 Rf_error(_("'%s' failed"), "cs_lusol");
94 UNPROTECT(1);
95 return b;
96}
97
98/* called from package MatrixModels's R code : */
99/* TODO: support NCOL(b) > 1 */
100/* TODO: give result list(L, coef, Xty, resid) */
101SEXP dgCMatrix_qrsol(SEXP a, SEXP b, SEXP order)
102{
103 /* FIXME? 'cs_qrsol' supports underdetermined systems. */
104 /* We should only require LENGTH(b) = max(m, n). */
105 int order_ = Rf_asInteger(order);
106 if (order_ < 0 || order_ > 3)
107 order_ = 0;
108 Matrix_cs *A = M2CXS(a, 1);
109 CXSPARSE_XTYPE_SET(CXSPARSE_REAL);
110 PROTECT(b = (TYPEOF(b) == REALSXP)
111 ? Rf_duplicate(b) : Rf_coerceVector(b, REALSXP));
112 if (LENGTH(b) != A->m)
113 Rf_error(_("dimensions of '%s' and '%s' are inconsistent"), "a", "b");
114 if (A->n <= 0 || A->n > A->m)
115 Rf_error(_("%s(%s, %s) requires m-by-n '%s' with m >= n > 0"),
116 "dgCMatrix_qrsol", "a", "b", "a");
117 if (!Matrix_cs_qrsol(order_, A, REAL(b)))
118 Rf_error(_("'%s' failed"), "cs_qrsol");
119 if (A->n < A->m) {
120 SEXP tmp = Rf_allocVector(REALSXP, A->n);
121 memcpy(REAL(tmp), REAL(b), sizeof(double) * (size_t) A->n);
122 b = tmp;
123 }
124 UNPROTECT(1);
125 return b;
126}
127
128/* called from package MatrixModels's R code : */
129/* TODO: support NCOL(b) > 1 */
130SEXP dgCMatrix_cholsol(SEXP at, SEXP b)
131{
132 /* Find least squares solution of A * X = B, given A' and B : */
133 cholmod_sparse *At = M2CHS(at, 1);
134 PROTECT(b = Rf_coerceVector(b, REALSXP));
135 if (LENGTH(b) != At->ncol)
136 Rf_error(_("dimensions of '%s' and '%s' are inconsistent"), "at", "b");
137 if (At->ncol <= 0 || At->ncol < At->nrow)
138 Rf_error(_("%s(%s, %s) requires m-by-n '%s' with n >= m > 0"),
139 "dgCMatrix_cholsol", "at", "b", "at");
140 double zero[] = { 0.0, 0.0 }, one[] = {1.0, 0.0}, mone[] = { -1.0, 0.0 };
141
142 /* L * L' = A' * A */
143 cholmod_factor *L = cholmod_analyze(At, &c);
144 if (!cholmod_factorize(At, L, &c))
145 Rf_error(_("'%s' failed"), "cholmod_factorize");
146
147 cholmod_dense *B = (cholmod_dense *) R_alloc(1, sizeof(cholmod_dense));
148 memset(B, 0, sizeof(cholmod_dense));
149 B->nrow = B->d = B->nzmax = (size_t) LENGTH(b);
150 B->ncol = 1;
151 B->x = REAL(b);
152 B->dtype = CHOLMOD_DOUBLE;
153 B->xtype = CHOLMOD_REAL;
154
155 /* A' * B = 1 * A' * B + 0 * <in/out> */
156 cholmod_dense *AtB = cholmod_allocate_dense(
157 At->nrow, 1, At->nrow, CHOLMOD_REAL, &c);
158 if (!cholmod_sdmult(At, 0, one, zero, B, AtB, &c))
159 Rf_error(_("'%s' failed"), "cholmod_sdmult");
160
161 /* C := solve(A' * A, A' * B) = solve(L', solve(L, A' * B)) */
162 cholmod_dense *C = cholmod_solve(CHOLMOD_A, L, AtB, &c);
163 if (!C)
164 Rf_error(_("'%s' failed"), "cholmod_solve");
165
166 /* R := A * A' * C - B = 1 * (A')' * A' * X + (-1) * B */
167 cholmod_dense *R = cholmod_copy_dense(B, &c);
168 if (!cholmod_sdmult(At, 1, mone, one, C, R, &c))
169 Rf_error(_("'%s' failed"), "cholmod_sdmult");
170
171 const char *nms[] = {"L", "coef", "Xty", "resid", ""};
172 SEXP ans = PROTECT(Rf_mkNamed(VECSXP, nms)), tmp;
173 /* L : */
174 PROTECT(tmp = CHF2M(L, 1));
175 SET_VECTOR_ELT(ans, 0, tmp);
176 /* coef : */
177 PROTECT(tmp = Rf_allocVector(REALSXP, (R_xlen_t) At->nrow));
178 memcpy(REAL(tmp), C->x, sizeof(double) * At->nrow);
179 SET_VECTOR_ELT(ans, 1, tmp);
180 /* Xty : */
181 PROTECT(tmp = Rf_allocVector(REALSXP, (R_xlen_t) At->nrow));
182 memcpy(REAL(tmp), AtB->x, sizeof(double) * At->nrow);
183 SET_VECTOR_ELT(ans, 2, tmp);
184 /* resid : */
185 PROTECT(tmp = Rf_allocVector(REALSXP, (R_xlen_t) At->ncol));
186 memcpy(REAL(tmp), R->x, sizeof(double) * At->ncol);
187 SET_VECTOR_ELT(ans, 3, tmp);
188
189 cholmod_free_factor(& L, &c);
190 cholmod_free_dense (&AtB, &c);
191 cholmod_free_dense (& C, &c);
192 cholmod_free_dense (& R, &c);
193
194 UNPROTECT(6);
195 return ans;
196}
197
198static
199int strmatch(const char *x, const char **valid)
200{
201 int i = 0;
202 while (valid[i][0] != '\0') {
203 if (strcmp(x, valid[i]) == 0)
204 return i;
205 ++i;
206 }
207 return -1;
208}
209
210/* <op>(diag(obj)) where obj=dCHMsimpl (LDLt) or obj=dtCMatrix (nonunit) */
211SEXP dtCMatrix_diag(SEXP obj, SEXP op)
212{
213 static const char *valid[] = {
214 /* 0 */ "trace",
215 /* 1 */ "sumLog",
216 /* 2 */ "prod",
217 /* 3 */ "min",
218 /* 4 */ "max",
219 /* 5 */ "range",
220 /* 6 */ "diag",
221 /* 7 */ "diagBack", ""};
222 int ivalid = -1;
223 if (TYPEOF(op) != STRSXP || LENGTH(op) < 1 ||
224 (op = STRING_ELT(op, 0)) == NA_STRING ||
225 (ivalid = strmatch(CHAR(op), valid)) < 0)
226 Rf_error(_("invalid '%s' to '%s'"), "op", __func__);
227
228 SEXP uplo = Rf_getAttrib(obj, Matrix_uploSym);
229 char ul = (TYPEOF(uplo) == STRSXP && LENGTH(uplo) > 0)
230 ? CHAR(STRING_ELT(uplo, 0))[0] : 'L';
231
232 SEXP p = PROTECT(GET_SLOT(obj, Matrix_pSym));
233 int *pp = INTEGER(p) + 1, j, k, kend, n = (int) (XLENGTH(p) - 1),
234 len = (ivalid < 5) ? 1 : ((ivalid < 6) ? 2 : n);
235
236 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym));
237 double *px = REAL(x), tmp;
238
239 SEXP ans = PROTECT(Rf_allocVector(REALSXP, len));
240 double *pans = REAL(ans);
241
242 switch (ivalid) {
243 case 0: /* trace */
244 pans[0] = 0.0;
245 for (j = 0, k = 0; j < n; ++j) {
246 kend = pp[j];
247 if (k < kend)
248 pans[0] += px[(ul == 'U') ? kend - 1 : k];
249 k = kend;
250 }
251 break;
252 case 1: /* sumLog */
253 pans[0] = 0.0;
254 for (j = 0, k = 0; j < n; ++j) {
255 kend = pp[j];
256 if (k < kend)
257 pans[0] += log(px[(ul == 'U') ? kend - 1 : k]);
258 k = kend;
259 }
260 break;
261 case 2: /* prod */
262 pans[0] = 1.0;
263 for (j = 0, k = 0; j < n; ++j) {
264 kend = pp[j];
265 if (k < kend)
266 pans[0] *= px[(ul == 'U') ? kend - 1 : k];
267 else
268 pans[0] *= 0.0;
269 k = kend;
270 }
271 break;
272 case 3: /* min */
273 pans[0] = R_PosInf;
274 for (j = 0, k = 0; j < n; ++j) {
275 kend = pp[j];
276 tmp = (k < kend) ? px[(ul == 'U') ? kend - 1 : k] : 0.0;
277 if (ISNAN(tmp)) {
278 pans[0] = tmp;
279 break;
280 }
281 else if (tmp < pans[0])
282 pans[0] = tmp;
283 k = kend;
284 }
285 break;
286 case 4: /* max */
287 pans[0] = R_NegInf;
288 for (j = 0, k = 0; j < n; ++j) {
289 kend = pp[j];
290 tmp = (k < kend) ? px[(ul == 'U') ? kend - 1 : k] : 0.0;
291 if (ISNAN(tmp)) {
292 pans[0] = tmp;
293 break;
294 }
295 else if (tmp > pans[0])
296 pans[0] = tmp;
297 k = kend;
298 }
299 break;
300 case 5: /* range */
301 pans[0] = R_PosInf;
302 pans[1] = R_NegInf;
303 for (j = 0, k = 0; j < n; ++j) {
304 kend = pp[j];
305 tmp = (k < kend) ? px[(ul == 'U') ? kend - 1 : k] : 0.0;
306 if (ISNAN(tmp)) {
307 pans[0] = pans[1] = tmp;
308 break;
309 }
310 else if (tmp < pans[0])
311 pans[0] = tmp;
312 else if (tmp > pans[1])
313 pans[1] = tmp;
314 k = kend;
315 }
316 break;
317 case 6: /* diag */
318 case 7: /* diagBack */
319 {
320
321 int *pperm = NULL;
322 if (ivalid == 7) {
323 SEXP perm = Rf_getAttrib(obj, Matrix_permSym);
324 if (TYPEOF(perm) == INTSXP && LENGTH(perm) == n)
325 pperm = INTEGER(perm);
326 }
327 for (j = 0, k = 0; j < n; ++j) {
328 kend = pp[j];
329 pans[(pperm) ? pperm[j] : j] =
330 (k < kend) ? px[(ul == 'U') ? kend - 1 : k] : 0.0;
331 k = kend;
332 }
333 break;
334 }
335 default:
336 break;
337 }
338
339 UNPROTECT(3);
340 return ans;
341}
342
343/* dmperm(x, nAns, seed) */
344SEXP Csparse_dmperm(SEXP x, SEXP nans, SEXP seed)
345{
346 Matrix_cs *A = M2CXS(x, 0);
347 CXSPARSE_XTYPE_SET(A->xtype);
348 Matrix_csd *D = Matrix_cs_dmperm(A, Rf_asInteger(seed));
349 if (!D)
350 return R_NilValue; /* MJ: why not an error ... ? */
351 int len = Rf_asInteger(nans);
352 if (len < 0)
353 len = 0;
354 else if (len > 6)
355 len = 6;
356 SEXP nms = PROTECT(Rf_allocVector(STRSXP, len)),
357 ans = PROTECT(Rf_allocVector(VECSXP, len)), tmp;
358 int k = len - 1;
359 switch (len) {
360 case 6:
361 SET_STRING_ELT(nms, k, Rf_mkChar("cc5"));
362 tmp = Rf_allocVector(INTSXP, 5);
363 memcpy(INTEGER(tmp), D->cc, sizeof(int) * 5);
364 SET_VECTOR_ELT(ans, k, tmp);
365 k--;
366 case 5:
367 SET_STRING_ELT(nms, k, Rf_mkChar("rr5"));
368 tmp = Rf_allocVector(INTSXP, 5);
369 memcpy(INTEGER(tmp), D->rr, sizeof(int) * 5);
370 SET_VECTOR_ELT(ans, k, tmp);
371 k--;
372 case 4:
373 SET_STRING_ELT(nms, k, Rf_mkChar("s"));
374 tmp = Rf_allocVector(INTSXP, D->nb + 1);
375 memcpy(INTEGER(tmp), D->s , sizeof(int) * (size_t) (D->nb + 1));
376 SET_VECTOR_ELT(ans, k, tmp);
377 k--;
378 case 3:
379 SET_STRING_ELT(nms, k, Rf_mkChar("r"));
380 tmp = Rf_allocVector(INTSXP, D->nb + 1);
381 memcpy(INTEGER(tmp), D->r , sizeof(int) * (size_t) (D->nb + 1));
382 SET_VECTOR_ELT(ans, k, tmp);
383 k--;
384 case 2:
385 SET_STRING_ELT(nms, k, Rf_mkChar("q"));
386 tmp = Rf_allocVector(INTSXP, A->n);
387 for (int j = 0, *q0 = D->q, *q1 = INTEGER(tmp); j < A->n; ++j)
388 q1[j] = q0[j] + 1;
389 SET_VECTOR_ELT(ans, k, tmp);
390 k--;
391 case 1:
392 SET_STRING_ELT(nms, k, Rf_mkChar("p"));
393 tmp = Rf_allocVector(INTSXP, A->m);
394 for (int i = 0, *p0 = D->p, *p1 = INTEGER(tmp); i < A->m; ++i)
395 p1[i] = p0[i] + 1;
396 SET_VECTOR_ELT(ans, k, tmp);
397 k--;
398 default:
399 break;
400 }
401 D = Matrix_cs_dfree(D);
402 Rf_setAttrib(ans, R_NamesSymbol, nms);
403 UNPROTECT(2);
404 return ans;
405}
406
407/* writeMM(obj, file) */
408SEXP Csparse_writeMM(SEXP obj, SEXP file)
409{
410 const char *class = Matrix_class(obj, valid_sparse_compressed, 6, __func__);
411
412 PROTECT_INDEX pid;
413 PROTECT_WITH_INDEX(obj, &pid);
414 if (class[0] == 'l' || class[1] == 'i') {
415 /* defined in ./coerce.c : */
416 SEXP sparse_as_kind(SEXP, const char *, char);
417 REPROTECT(obj = sparse_as_kind(obj, class, 'd'), pid);
418 class = Matrix_class(obj, valid_sparse_compressed, 6, __func__);
419 }
420 if (class[1] == 't') {
421 /* defined in ./coerce.c : */
422 SEXP sparse_as_general(SEXP, const char *);
423 REPROTECT(obj = sparse_as_general(obj, class), pid);
424 class = Matrix_class(obj, valid_sparse_compressed, 6, __func__);
425 }
426
427 cholmod_sparse *A = M2CHS(obj, 1);
428 const char *filename = CHAR(Rf_asChar(file));
429 FILE *f = fopen(filename, "w");
430 if (!f)
431 Rf_error(_("failed to open file \"%s\" for writing"), filename);
432 if (!cholmod_write_sparse(f, A, (cholmod_sparse *) NULL, (char *) NULL, &c))
433 Rf_error(_("'%s' failed"), "cholmod_write_sparse");
434 fclose(f);
435
436 UNPROTECT(1);
437 return R_NilValue;
438}
sparse_as_general
SEXP sparse_as_general(SEXP, const char *)
Definition coerce.c:2298
sparse_as_kind
SEXP sparse_as_kind(SEXP, const char *, char)
Definition coerce.c:2076
strmatch
static int strmatch(const char *x, const char **valid)
Definition Csparse.c:199
dtCMatrix_diag
SEXP dtCMatrix_diag(SEXP obj, SEXP op)
Definition Csparse.c:211
Csparse_dmperm
SEXP Csparse_dmperm(SEXP x, SEXP nans, SEXP seed)
Definition Csparse.c:344
checkpi
SEXP checkpi(SEXP dim, SEXP p, SEXP i)
Definition Csparse.c:6
dgCMatrix_qrsol
SEXP dgCMatrix_qrsol(SEXP a, SEXP b, SEXP order)
Definition Csparse.c:101
dgCMatrix_lusol
SEXP dgCMatrix_lusol(SEXP a, SEXP b)
Definition Csparse.c:82
Csparse_writeMM
SEXP Csparse_writeMM(SEXP obj, SEXP file)
Definition Csparse.c:408
CsparseMatrix_validate_maybe_sorting
SEXP CsparseMatrix_validate_maybe_sorting(SEXP x)
Definition Csparse.c:64
dgCMatrix_cholsol
SEXP dgCMatrix_cholsol(SEXP at, SEXP b)
Definition Csparse.c:130
_
#define _(String)
Definition Mdefines.h:66
Matrix_class
const char * Matrix_class(SEXP, const char **, int, const char *)
Definition objects.c:112
GET_SLOT
#define GET_SLOT(x, name)
Definition Mdefines.h:72
errorString
#define errorString(...)
Definition Mdefines.h:69
valid_sparse_compressed
const char * valid_sparse_compressed[]
Definition Mdefines.h:329
TYPEOF
#define TYPEOF(s)
Definition Mdefines.h:123
M2CHS
cholmod_sparse * M2CHS(SEXP obj, int values)
Definition cholmod-etc.c:132
CHF2M
SEXP CHF2M(cholmod_factor *L, int values)
Definition cholmod-etc.c:222
c
cholmod_common c
Definition cholmod-etc.c:5
Matrix_cs_qrsol
int Matrix_cs_qrsol(int order, const Matrix_cs *A, void *b)
Definition cs-etc.c:258
Matrix_cs_dmperm
Matrix_csd * Matrix_cs_dmperm(const Matrix_cs *A, int seed)
Definition cs-etc.c:89
Matrix_cs_lusol
int Matrix_cs_lusol(int order, const Matrix_cs *A, void *b, double tol)
Definition cs-etc.c:173
Matrix_cs_dfree
Matrix_csd * Matrix_cs_dfree(Matrix_csd *D)
Definition cs-etc.c:77
M2CXS
Matrix_cs * M2CXS(SEXP obj, int values)
Definition cs-etc.c:6
CXSPARSE_XTYPE_SET
#define CXSPARSE_XTYPE_SET(_VALUE_)
Definition cs-etc.h:12
CXSPARSE_REAL
#define CXSPARSE_REAL
Definition cs-etc.h:8
Matrix_permSym
SEXP Matrix_permSym
Definition init.c:623
Matrix_DimSym
SEXP Matrix_DimSym
Definition init.c:598
Matrix_xSym
SEXP Matrix_xSym
Definition init.c:635
Matrix_iSym
SEXP Matrix_iSym
Definition init.c:607
Matrix_uploSym
SEXP Matrix_uploSym
Definition init.c:632
Matrix_pSym
SEXP Matrix_pSym
Definition init.c:622
Matrix_cs_dmperm_results
Definition cs-etc.h:49
Matrix_cs_dmperm_results::rr
int rr[5]
Definition cs-etc.h:55
Matrix_cs_dmperm_results::r
int * r
Definition cs-etc.h:52
Matrix_cs_dmperm_results::s
int * s
Definition cs-etc.h:53
Matrix_cs_dmperm_results::nb
int nb
Definition cs-etc.h:54
Matrix_cs_dmperm_results::q
int * q
Definition cs-etc.h:51
Matrix_cs_dmperm_results::p
int * p
Definition cs-etc.h:50
Matrix_cs_dmperm_results::cc
int cc[5]
Definition cs-etc.h:56
Matrix_cs_sparse
Definition cs-etc.h:17
Matrix_cs_sparse::n
int n
Definition cs-etc.h:20
Matrix_cs_sparse::m
int m
Definition cs-etc.h:19
Matrix_cs_sparse::xtype
int xtype
Definition cs-etc.h:25
Generated by doxygen 1.12.0