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