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