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