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