SCM Repository
Annotation of /pkg/src/Csparse.c
Parent Directory
|
Revision Log
Revision 1660 - (view) (download) (as text)
1 : | bates | 1218 | /* Sparse matrices in compressed column-oriented form */ |
2 : | bates | 922 | #include "Csparse.h" |
3 : | #include "chm_common.h" | ||
4 : | |||
5 : | SEXP Csparse_validate(SEXP x) | ||
6 : | { | ||
7 : | maechler | 1575 | /* NB: we do *NOT* check a potential 'x' slot here, at all */ |
8 : | bates | 922 | SEXP pslot = GET_SLOT(x, Matrix_pSym), |
9 : | islot = GET_SLOT(x, Matrix_iSym); | ||
10 : | maechler | 1660 | int j, k, sorted, |
11 : | bates | 922 | *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), |
12 : | maechler | 1660 | nrow = dims[0], |
13 : | ncol = dims[1], | ||
14 : | maechler | 1654 | *xp = INTEGER(pslot), |
15 : | bates | 922 | *xi = INTEGER(islot); |
16 : | |||
17 : | maechler | 1654 | if (length(pslot) != dims[1] + 1) |
18 : | return mkString(_("slot p must have length = ncol(.) + 1")); | ||
19 : | bates | 922 | if (xp[0] != 0) |
20 : | return mkString(_("first element of slot p must be zero")); | ||
21 : | if (length(islot) != xp[ncol]) | ||
22 : | bates | 1555 | return |
23 : | mkString(_("last element of slot p must match length of slots i and x")); | ||
24 : | for (j = 0; j < length(islot); j++) { | ||
25 : | if (xi[j] < 0 || xi[j] >= nrow) | ||
26 : | return mkString(_("all row indices must be between 0 and nrow-1")); | ||
27 : | } | ||
28 : | sorted = TRUE; | ||
29 : | bates | 922 | for (j = 0; j < ncol; j++) { |
30 : | if (xp[j] > xp[j+1]) | ||
31 : | return mkString(_("slot p must be non-decreasing")); | ||
32 : | bates | 1555 | for (k = xp[j] + 1; k < xp[j + 1]; k++) |
33 : | if (xi[k] < xi[k - 1]) sorted = FALSE; | ||
34 : | bates | 922 | } |
35 : | maechler | 1654 | if (!sorted) { |
36 : | cholmod_sparse *chx = as_cholmod_sparse(x); | ||
37 : | cholmod_sort(chx, &c); | ||
38 : | Free(chx); | ||
39 : | } | ||
40 : | bates | 922 | return ScalarLogical(1); |
41 : | } | ||
42 : | |||
43 : | bates | 1059 | SEXP Csparse_to_dense(SEXP x) |
44 : | { | ||
45 : | cholmod_sparse *chxs = as_cholmod_sparse(x); | ||
46 : | cholmod_dense *chxd = cholmod_sparse_to_dense(chxs, &c); | ||
47 : | |||
48 : | bates | 1141 | Free(chxs); |
49 : | maechler | 1660 | return chm_dense_to_SEXP(chxd, 1, Real_kind(x), |
50 : | GET_SLOT(x, Matrix_DimNamesSym)); | ||
51 : | bates | 1059 | } |
52 : | |||
53 : | maechler | 1548 | SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri) |
54 : | bates | 1371 | { |
55 : | cholmod_sparse *chxs = as_cholmod_sparse(x); | ||
56 : | cholmod_sparse | ||
57 : | *chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c); | ||
58 : | int uploT = 0; char *diag = ""; | ||
59 : | |||
60 : | Free(chxs); | ||
61 : | if (asLogical(tri)) { /* triangular sparse matrices */ | ||
62 : | uploT = (strcmp(CHAR(asChar(GET_SLOT(x, Matrix_uploSym))), "U")) ? | ||
63 : | -1 : 1; | ||
64 : | diag = CHAR(asChar(GET_SLOT(x, Matrix_diagSym))); | ||
65 : | } | ||
66 : | maechler | 1548 | return chm_sparse_to_SEXP(chxcp, 1, uploT, 0, diag, |
67 : | GET_SLOT(x, Matrix_DimNamesSym)); | ||
68 : | bates | 1371 | } |
69 : | |||
70 : | bates | 1366 | SEXP Csparse_to_matrix(SEXP x) |
71 : | bates | 922 | { |
72 : | maechler | 925 | cholmod_sparse *chxs = as_cholmod_sparse(x); |
73 : | bates | 1366 | cholmod_dense *chxd = cholmod_sparse_to_dense(chxs, &c); |
74 : | |||
75 : | Free(chxs); | ||
76 : | return chm_dense_to_matrix(chxd, 1, | ||
77 : | GET_SLOT(x, Matrix_DimNamesSym)); | ||
78 : | } | ||
79 : | |||
80 : | SEXP Csparse_to_Tsparse(SEXP x, SEXP tri) | ||
81 : | { | ||
82 : | cholmod_sparse *chxs = as_cholmod_sparse(x); | ||
83 : | bates | 922 | cholmod_triplet *chxt = cholmod_sparse_to_triplet(chxs, &c); |
84 : | maechler | 1548 | int uploT = 0; |
85 : | char *diag = ""; | ||
86 : | int Rkind = (chxs->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; | ||
87 : | bates | 922 | |
88 : | bates | 1141 | Free(chxs); |
89 : | bates | 1366 | if (asLogical(tri)) { /* triangular sparse matrices */ |
90 : | bates | 1419 | uploT = (*uplo_P(x) == 'U') ? -1 : 1; |
91 : | diag = diag_P(x); | ||
92 : | bates | 1366 | } |
93 : | maechler | 1548 | return chm_triplet_to_SEXP(chxt, 1, uploT, Rkind, diag, |
94 : | bates | 1371 | GET_SLOT(x, Matrix_DimNamesSym)); |
95 : | bates | 922 | } |
96 : | |||
97 : | bates | 1448 | /* this used to be called sCMatrix_to_gCMatrix(..) [in ./dsCMatrix.c ]: */ |
98 : | bates | 1371 | SEXP Csparse_symmetric_to_general(SEXP x) |
99 : | { | ||
100 : | cholmod_sparse *chx = as_cholmod_sparse(x), *chgx; | ||
101 : | maechler | 1548 | int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
102 : | bates | 1371 | |
103 : | if (!(chx->stype)) | ||
104 : | maechler | 1548 | error(_("Nonsymmetric matrix in Csparse_symmetric_to_general")); |
105 : | maechler | 1375 | chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c); |
106 : | /* xtype: pattern, "real", complex or .. */ | ||
107 : | bates | 1371 | Free(chx); |
108 : | maechler | 1548 | return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
109 : | bates | 1371 | GET_SLOT(x, Matrix_DimNamesSym)); |
110 : | } | ||
111 : | |||
112 : | maechler | 1618 | SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo) |
113 : | maechler | 1598 | { |
114 : | cholmod_sparse *chx = as_cholmod_sparse(x), *chgx; | ||
115 : | maechler | 1618 | int uploT = (*CHAR(asChar(uplo)) == 'U') ? -1 : 1; |
116 : | maechler | 1598 | int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
117 : | |||
118 : | maechler | 1618 | chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c); |
119 : | maechler | 1598 | /* xtype: pattern, "real", complex or .. */ |
120 : | Free(chx); | ||
121 : | return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", | ||
122 : | GET_SLOT(x, Matrix_DimNamesSym)); | ||
123 : | } | ||
124 : | |||
125 : | bates | 1369 | SEXP Csparse_transpose(SEXP x, SEXP tri) |
126 : | bates | 922 | { |
127 : | cholmod_sparse *chx = as_cholmod_sparse(x); | ||
128 : | maechler | 1568 | int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
129 : | bates | 922 | cholmod_sparse *chxt = cholmod_transpose(chx, (int) chx->xtype, &c); |
130 : | bates | 1366 | SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp; |
131 : | bates | 1369 | int uploT = 0; char *diag = ""; |
132 : | |||
133 : | bates | 1141 | Free(chx); |
134 : | bates | 1366 | tmp = VECTOR_ELT(dn, 0); /* swap the dimnames */ |
135 : | SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1)); | ||
136 : | SET_VECTOR_ELT(dn, 1, tmp); | ||
137 : | UNPROTECT(1); | ||
138 : | bates | 1369 | if (asLogical(tri)) { /* triangular sparse matrices */ |
139 : | bates | 1419 | uploT = (*uplo_P(x) == 'U') ? -1 : 1; |
140 : | diag = diag_P(x); | ||
141 : | bates | 1369 | } |
142 : | maechler | 1548 | return chm_sparse_to_SEXP(chxt, 1, uploT, Rkind, diag, dn); |
143 : | bates | 922 | } |
144 : | |||
145 : | SEXP Csparse_Csparse_prod(SEXP a, SEXP b) | ||
146 : | { | ||
147 : | maechler | 1659 | cholmod_sparse |
148 : | *cha = as_cholmod_sparse(a), | ||
149 : | bates | 1059 | *chb = as_cholmod_sparse(b); |
150 : | cholmod_sparse *chc = cholmod_ssmult(cha, chb, 0, cha->xtype, 1, &c); | ||
151 : | bates | 1366 | SEXP dn = allocVector(VECSXP, 2); |
152 : | bates | 922 | |
153 : | bates | 1141 | Free(cha); Free(chb); |
154 : | bates | 1366 | SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
155 : | duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); | ||
156 : | SET_VECTOR_ELT(dn, 1, | ||
157 : | duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1))); | ||
158 : | maechler | 1548 | return chm_sparse_to_SEXP(chc, 1, 0, 0, "", dn); |
159 : | bates | 922 | } |
160 : | |||
161 : | maechler | 1659 | SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b, SEXP trans) |
162 : | bates | 1657 | { |
163 : | maechler | 1659 | int tr = asLogical(trans); |
164 : | cholmod_sparse | ||
165 : | *cha = as_cholmod_sparse(a), | ||
166 : | bates | 1657 | *chb = as_cholmod_sparse(b); |
167 : | maechler | 1659 | cholmod_sparse *chTr, *chc; |
168 : | bates | 1657 | SEXP dn = allocVector(VECSXP, 2); |
169 : | |||
170 : | maechler | 1659 | /* cholmod_sparse *chTr = cholmod_transpose(cha, 1, &c); */ |
171 : | /* cholmod_sparse *chc = cholmod_ssmult(chTr, chb, 0, cha->xtype, 1, &c); */ | ||
172 : | |||
173 : | if (tr) | ||
174 : | chTr = cholmod_transpose(chb, chb->xtype, &c); | ||
175 : | else | ||
176 : | chTr = cholmod_transpose(cha, cha->xtype, &c); | ||
177 : | chc = cholmod_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb, | ||
178 : | 0, cha->xtype, 1, &c); | ||
179 : | |||
180 : | Free(cha); Free(chb); cholmod_free_sparse(&chTr, &c); | ||
181 : | |||
182 : | bates | 1657 | SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
183 : | maechler | 1659 | duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), (tr) ? 0 : 1))); |
184 : | bates | 1657 | SET_VECTOR_ELT(dn, 1, |
185 : | maechler | 1659 | duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), (tr) ? 0 : 1))); |
186 : | bates | 1657 | return chm_sparse_to_SEXP(chc, 1, 0, 0, "", dn); |
187 : | } | ||
188 : | |||
189 : | bates | 922 | SEXP Csparse_dense_prod(SEXP a, SEXP b) |
190 : | { | ||
191 : | cholmod_sparse *cha = as_cholmod_sparse(a); | ||
192 : | maechler | 1660 | SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
193 : | cholmod_dense *chb = as_cholmod_dense(b_M); | ||
194 : | maechler | 1548 | cholmod_dense *chc = |
195 : | bates | 1448 | cholmod_allocate_dense(cha->nrow, chb->ncol, cha->nrow, chb->xtype, &c); |
196 : | maechler | 1660 | SEXP dn = allocVector(VECSXP, 2); |
197 : | bates | 1448 | double alpha[] = {1,0}, beta[] = {0,0}; |
198 : | bates | 922 | |
199 : | bates | 1448 | cholmod_sdmult(cha, 0, alpha, beta, chb, chc, &c); |
200 : | bates | 1141 | Free(cha); Free(chb); |
201 : | bates | 1461 | UNPROTECT(1); |
202 : | maechler | 1660 | SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
203 : | duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); | ||
204 : | SET_VECTOR_ELT(dn, 1, | ||
205 : | duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1))); | ||
206 : | return chm_dense_to_SEXP(chc, 1, 0, dn); | ||
207 : | bates | 922 | } |
208 : | maechler | 925 | |
209 : | bates | 1067 | SEXP Csparse_dense_crossprod(SEXP a, SEXP b) |
210 : | { | ||
211 : | cholmod_sparse *cha = as_cholmod_sparse(a); | ||
212 : | maechler | 1660 | SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
213 : | cholmod_dense *chb = as_cholmod_dense(b_M); | ||
214 : | bates | 1448 | cholmod_dense *chc = |
215 : | cholmod_allocate_dense(cha->ncol, chb->ncol, cha->ncol, chb->xtype, &c); | ||
216 : | maechler | 1660 | SEXP dn = allocVector(VECSXP, 2); |
217 : | bates | 1448 | double alpha[] = {1,0}, beta[] = {0,0}; |
218 : | bates | 1067 | |
219 : | bates | 1448 | cholmod_sdmult(cha, 1, alpha, beta, chb, chc, &c); |
220 : | bates | 1067 | Free(cha); Free(chb); |
221 : | bates | 1461 | UNPROTECT(1); |
222 : | maechler | 1660 | SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
223 : | duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1))); | ||
224 : | SET_VECTOR_ELT(dn, 1, | ||
225 : | duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1))); | ||
226 : | return chm_dense_to_SEXP(chc, 1, 0, dn); | ||
227 : | bates | 1067 | } |
228 : | |||
229 : | maechler | 1659 | /* Computes x'x or x x' -- see Csparse_Csparse_crossprod above for x'y and x y' */ |
230 : | bates | 928 | SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet) |
231 : | bates | 922 | { |
232 : | maechler | 957 | int trip = asLogical(triplet), |
233 : | tr = asLogical(trans); /* gets reversed because _aat is tcrossprod */ | ||
234 : | bates | 928 | cholmod_triplet |
235 : | *cht = trip ? as_cholmod_triplet(x) : (cholmod_triplet*) NULL; | ||
236 : | cholmod_sparse *chcp, *chxt, | ||
237 : | *chx = trip ? cholmod_triplet_to_sparse(cht, cht->nnz, &c) | ||
238 : | : as_cholmod_sparse(x); | ||
239 : | bates | 1366 | SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
240 : | bates | 922 | |
241 : | bates | 923 | if (!tr) |
242 : | bates | 1353 | chxt = cholmod_transpose(chx, chx->xtype, &c); |
243 : | bates | 928 | chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c); |
244 : | maechler | 957 | if(!chcp) |
245 : | maechler | 1618 | error(_("Csparse_crossprod(): error return from cholmod_aat()")); |
246 : | bates | 1360 | cholmod_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c); |
247 : | chcp->stype = 1; | ||
248 : | bates | 930 | if (trip) { |
249 : | cholmod_free_sparse(&chx, &c); | ||
250 : | bates | 1141 | Free(cht); |
251 : | bates | 930 | } else { |
252 : | bates | 1141 | Free(chx); |
253 : | bates | 930 | } |
254 : | bates | 923 | if (!tr) cholmod_free_sparse(&chxt, &c); |
255 : | bates | 1366 | /* create dimnames */ |
256 : | SET_VECTOR_ELT(dn, 0, | ||
257 : | duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), | ||
258 : | maechler | 1660 | (tr) ? 0 : 1))); |
259 : | bates | 1366 | SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(dn, 0))); |
260 : | UNPROTECT(1); | ||
261 : | maechler | 1548 | return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn); |
262 : | bates | 922 | } |
263 : | bates | 923 | |
264 : | maechler | 1618 | SEXP Csparse_drop(SEXP x, SEXP tol) |
265 : | { | ||
266 : | cholmod_sparse *chx = as_cholmod_sparse(x), | ||
267 : | *ans = cholmod_copy(chx, chx->stype, chx->xtype, &c); | ||
268 : | double dtol = asReal(tol); | ||
269 : | int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; | ||
270 : | |||
271 : | if(!cholmod_drop(dtol, ans, &c)) | ||
272 : | error(_("cholmod_drop() failed")); | ||
273 : | Free(chx); | ||
274 : | /* FIXME: currently drops dimnames */ | ||
275 : | return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue); | ||
276 : | } | ||
277 : | |||
278 : | |||
279 : | bates | 1218 | SEXP Csparse_horzcat(SEXP x, SEXP y) |
280 : | { | ||
281 : | cholmod_sparse *chx = as_cholmod_sparse(x), | ||
282 : | *chy = as_cholmod_sparse(y), *ans; | ||
283 : | maechler | 1548 | int Rkind = 0; /* only for "d" - FIXME */ |
284 : | maechler | 1375 | |
285 : | bates | 1218 | ans = cholmod_horzcat(chx, chy, 1, &c); |
286 : | Free(chx); Free(chy); | ||
287 : | bates | 1366 | /* FIXME: currently drops dimnames */ |
288 : | maechler | 1548 | return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue); |
289 : | bates | 1218 | } |
290 : | |||
291 : | SEXP Csparse_vertcat(SEXP x, SEXP y) | ||
292 : | { | ||
293 : | cholmod_sparse *chx = as_cholmod_sparse(x), | ||
294 : | *chy = as_cholmod_sparse(y), *ans; | ||
295 : | maechler | 1548 | int Rkind = 0; /* only for "d" - FIXME */ |
296 : | maechler | 1375 | |
297 : | bates | 1218 | ans = cholmod_vertcat(chx, chy, 1, &c); |
298 : | Free(chx); Free(chy); | ||
299 : | bates | 1366 | /* FIXME: currently drops dimnames */ |
300 : | maechler | 1548 | return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue); |
301 : | bates | 1218 | } |
302 : | bates | 1265 | |
303 : | SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2) | ||
304 : | { | ||
305 : | cholmod_sparse *chx = as_cholmod_sparse(x), *ans; | ||
306 : | maechler | 1548 | int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
307 : | bates | 1265 | |
308 : | ans = cholmod_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c); | ||
309 : | Free(chx); | ||
310 : | maechler | 1548 | return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue); |
311 : | bates | 1265 | } |
312 : | bates | 1366 | |
313 : | SEXP Csparse_diagU2N(SEXP x) | ||
314 : | { | ||
315 : | cholmod_sparse *chx = as_cholmod_sparse(x); | ||
316 : | cholmod_sparse *eye = cholmod_speye(chx->nrow, chx->ncol, chx->xtype, &c); | ||
317 : | double one[] = {1, 0}; | ||
318 : | cholmod_sparse *ans = cholmod_add(chx, eye, one, one, TRUE, TRUE, &c); | ||
319 : | int uploT = (strcmp(CHAR(asChar(GET_SLOT(x, Matrix_uploSym))), "U")) ? | ||
320 : | -1 : 1; | ||
321 : | maechler | 1548 | int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
322 : | bates | 1366 | |
323 : | Free(chx); cholmod_free_sparse(&eye, &c); | ||
324 : | maechler | 1548 | return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N", |
325 : | bates | 1366 | duplicate(GET_SLOT(x, Matrix_DimNamesSym))); |
326 : | } | ||
327 : | |||
328 : | SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j) | ||
329 : | { | ||
330 : | cholmod_sparse *chx = as_cholmod_sparse(x); | ||
331 : | int rsize = (isNull(i)) ? -1 : LENGTH(i), | ||
332 : | csize = (isNull(j)) ? -1 : LENGTH(j); | ||
333 : | maechler | 1548 | int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
334 : | bates | 1366 | |
335 : | if (rsize >= 0 && !isInteger(i)) | ||
336 : | error(_("Index i must be NULL or integer")); | ||
337 : | if (csize >= 0 && !isInteger(j)) | ||
338 : | error(_("Index j must be NULL or integer")); | ||
339 : | return chm_sparse_to_SEXP(cholmod_submatrix(chx, INTEGER(i), rsize, | ||
340 : | maechler | 1375 | INTEGER(j), csize, |
341 : | bates | 1366 | TRUE, TRUE, &c), |
342 : | maechler | 1548 | 1, 0, Rkind, "", R_NilValue); |
343 : | bates | 1366 | } |
R-Forge@R-project.org | ViewVC Help |
Powered by ViewVC 1.0.0 |