SCM Repository
Annotation of /pkg/src/Csparse.c
Parent Directory
|
Revision Log
Revision 2120 - (view) (download) (as text)
1 : | bates | 1218 | /* Sparse matrices in compressed column-oriented form */ |
2 : | bates | 922 | #include "Csparse.h" |
3 : | maechler | 2120 | #include "Tsparse.h" |
4 : | bates | 922 | #include "chm_common.h" |
5 : | |||
6 : | SEXP Csparse_validate(SEXP x) | ||
7 : | { | ||
8 : | maechler | 1575 | /* NB: we do *NOT* check a potential 'x' slot here, at all */ |
9 : | bates | 922 | SEXP pslot = GET_SLOT(x, Matrix_pSym), |
10 : | islot = GET_SLOT(x, Matrix_iSym); | ||
11 : | maechler | 1893 | Rboolean sorted, strictly; |
12 : | int j, k, | ||
13 : | bates | 922 | *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), |
14 : | maechler | 1660 | nrow = dims[0], |
15 : | ncol = dims[1], | ||
16 : | maechler | 1654 | *xp = INTEGER(pslot), |
17 : | bates | 922 | *xi = INTEGER(islot); |
18 : | |||
19 : | maechler | 1654 | if (length(pslot) != dims[1] + 1) |
20 : | return mkString(_("slot p must have length = ncol(.) + 1")); | ||
21 : | bates | 922 | if (xp[0] != 0) |
22 : | return mkString(_("first element of slot p must be zero")); | ||
23 : | maechler | 1893 | if (length(islot) < xp[ncol]) /* allow larger slots from over-allocation!*/ |
24 : | bates | 1555 | return |
25 : | mkString(_("last element of slot p must match length of slots i and x")); | ||
26 : | for (j = 0; j < length(islot); j++) { | ||
27 : | if (xi[j] < 0 || xi[j] >= nrow) | ||
28 : | return mkString(_("all row indices must be between 0 and nrow-1")); | ||
29 : | } | ||
30 : | maechler | 1893 | sorted = TRUE; strictly = TRUE; |
31 : | bates | 922 | for (j = 0; j < ncol; j++) { |
32 : | if (xp[j] > xp[j+1]) | ||
33 : | return mkString(_("slot p must be non-decreasing")); | ||
34 : | maechler | 1893 | if(sorted) |
35 : | for (k = xp[j] + 1; k < xp[j + 1]; k++) { | ||
36 : | if (xi[k] < xi[k - 1]) | ||
37 : | sorted = FALSE; | ||
38 : | else if (xi[k] == xi[k - 1]) | ||
39 : | strictly = FALSE; | ||
40 : | } | ||
41 : | bates | 922 | } |
42 : | maechler | 1654 | if (!sorted) { |
43 : | maechler | 1960 | CHM_SP chx = AS_CHM_SP(x); |
44 : | R_CheckStack(); | ||
45 : | |||
46 : | maechler | 1654 | cholmod_sort(chx, &c); |
47 : | maechler | 1893 | /* Now re-check that row indices are *strictly* increasing |
48 : | * (and not just increasing) within each column : */ | ||
49 : | for (j = 0; j < ncol; j++) { | ||
50 : | for (k = xp[j] + 1; k < xp[j + 1]; k++) | ||
51 : | if (xi[k] == xi[k - 1]) | ||
52 : | return mkString(_("slot i is not *strictly* increasing inside a column (even after cholmod_sort)")); | ||
53 : | } | ||
54 : | |||
55 : | } else if(!strictly) { /* sorted, but not strictly */ | ||
56 : | return mkString(_("slot i is not *strictly* increasing inside a column")); | ||
57 : | maechler | 1654 | } |
58 : | bates | 922 | return ScalarLogical(1); |
59 : | } | ||
60 : | |||
61 : | maechler | 1968 | SEXP Rsparse_validate(SEXP x) |
62 : | { | ||
63 : | /* NB: we do *NOT* check a potential 'x' slot here, at all */ | ||
64 : | SEXP pslot = GET_SLOT(x, Matrix_pSym), | ||
65 : | jslot = GET_SLOT(x, Matrix_jSym); | ||
66 : | Rboolean sorted, strictly; | ||
67 : | int i, k, | ||
68 : | *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), | ||
69 : | nrow = dims[0], | ||
70 : | ncol = dims[1], | ||
71 : | *xp = INTEGER(pslot), | ||
72 : | *xj = INTEGER(jslot); | ||
73 : | |||
74 : | if (length(pslot) != dims[0] + 1) | ||
75 : | return mkString(_("slot p must have length = nrow(.) + 1")); | ||
76 : | if (xp[0] != 0) | ||
77 : | return mkString(_("first element of slot p must be zero")); | ||
78 : | if (length(jslot) < xp[nrow]) /* allow larger slots from over-allocation!*/ | ||
79 : | return | ||
80 : | mkString(_("last element of slot p must match length of slots j and x")); | ||
81 : | for (i = 0; i < length(jslot); i++) { | ||
82 : | if (xj[i] < 0 || xj[i] >= ncol) | ||
83 : | return mkString(_("all column indices must be between 0 and ncol-1")); | ||
84 : | } | ||
85 : | sorted = TRUE; strictly = TRUE; | ||
86 : | for (i = 0; i < nrow; i++) { | ||
87 : | if (xp[i] > xp[i+1]) | ||
88 : | return mkString(_("slot p must be non-decreasing")); | ||
89 : | if(sorted) | ||
90 : | for (k = xp[i] + 1; k < xp[i + 1]; k++) { | ||
91 : | if (xj[k] < xj[k - 1]) | ||
92 : | sorted = FALSE; | ||
93 : | else if (xj[k] == xj[k - 1]) | ||
94 : | strictly = FALSE; | ||
95 : | } | ||
96 : | } | ||
97 : | if (!sorted) | ||
98 : | /* cannot easily use cholmod_sort(.) ... -> "error out" :*/ | ||
99 : | return mkString(_("slot j is not increasing inside a column")); | ||
100 : | else if(!strictly) /* sorted, but not strictly */ | ||
101 : | return mkString(_("slot j is not *strictly* increasing inside a column")); | ||
102 : | |||
103 : | return ScalarLogical(1); | ||
104 : | } | ||
105 : | |||
106 : | |||
107 : | maechler | 1751 | /* Called from ../R/Csparse.R : */ |
108 : | /* Can only return [dln]geMatrix (no symm/triang); | ||
109 : | * FIXME: replace by non-CHOLMOD code ! */ | ||
110 : | bates | 1059 | SEXP Csparse_to_dense(SEXP x) |
111 : | { | ||
112 : | maechler | 1960 | CHM_SP chxs = AS_CHM_SP(x); |
113 : | maechler | 1751 | /* This loses the symmetry property, since cholmod_dense has none, |
114 : | * BUT, much worse (FIXME!), it also transforms CHOLMOD_PATTERN ("n") matrices | ||
115 : | * to numeric (CHOLMOD_REAL) ones : */ | ||
116 : | maechler | 1960 | CHM_DN chxd = cholmod_sparse_to_dense(chxs, &c); |
117 : | maechler | 1751 | int Rkind = (chxs->xtype == CHOLMOD_PATTERN)? -1 : Real_kind(x); |
118 : | maechler | 1960 | R_CheckStack(); |
119 : | bates | 1059 | |
120 : | maechler | 1736 | return chm_dense_to_SEXP(chxd, 1, Rkind, GET_SLOT(x, Matrix_DimNamesSym)); |
121 : | bates | 1059 | } |
122 : | |||
123 : | maechler | 1548 | SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri) |
124 : | bates | 1371 | { |
125 : | maechler | 1960 | CHM_SP chxs = AS_CHM_SP(x); |
126 : | CHM_SP chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c); | ||
127 : | bates | 1867 | int tr = asLogical(tri); |
128 : | maechler | 1960 | R_CheckStack(); |
129 : | bates | 1371 | |
130 : | maechler | 1960 | return chm_sparse_to_SEXP(chxcp, 1/*do_free*/, |
131 : | bates | 1867 | tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, |
132 : | 0, tr ? diag_P(x) : "", | ||
133 : | maechler | 1548 | GET_SLOT(x, Matrix_DimNamesSym)); |
134 : | bates | 1371 | } |
135 : | |||
136 : | bates | 1366 | SEXP Csparse_to_matrix(SEXP x) |
137 : | bates | 922 | { |
138 : | maechler | 1960 | return chm_dense_to_matrix(cholmod_sparse_to_dense(AS_CHM_SP(x), &c), |
139 : | 1 /*do_free*/, GET_SLOT(x, Matrix_DimNamesSym)); | ||
140 : | bates | 1366 | } |
141 : | |||
142 : | SEXP Csparse_to_Tsparse(SEXP x, SEXP tri) | ||
143 : | { | ||
144 : | maechler | 1960 | CHM_SP chxs = AS_CHM_SP(x); |
145 : | CHM_TR chxt = cholmod_sparse_to_triplet(chxs, &c); | ||
146 : | bates | 1867 | int tr = asLogical(tri); |
147 : | maechler | 1736 | int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
148 : | maechler | 1960 | R_CheckStack(); |
149 : | bates | 922 | |
150 : | bates | 1867 | return chm_triplet_to_SEXP(chxt, 1, |
151 : | tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, | ||
152 : | Rkind, tr ? diag_P(x) : "", | ||
153 : | bates | 1371 | GET_SLOT(x, Matrix_DimNamesSym)); |
154 : | bates | 922 | } |
155 : | |||
156 : | bates | 1448 | /* this used to be called sCMatrix_to_gCMatrix(..) [in ./dsCMatrix.c ]: */ |
157 : | bates | 1371 | SEXP Csparse_symmetric_to_general(SEXP x) |
158 : | { | ||
159 : | maechler | 1960 | CHM_SP chx = AS_CHM_SP(x), chgx; |
160 : | maechler | 1736 | int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
161 : | maechler | 1960 | R_CheckStack(); |
162 : | bates | 1371 | |
163 : | if (!(chx->stype)) | ||
164 : | maechler | 1548 | error(_("Nonsymmetric matrix in Csparse_symmetric_to_general")); |
165 : | maechler | 1375 | chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c); |
166 : | /* xtype: pattern, "real", complex or .. */ | ||
167 : | maechler | 1548 | return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
168 : | bates | 1371 | GET_SLOT(x, Matrix_DimNamesSym)); |
169 : | } | ||
170 : | |||
171 : | maechler | 1618 | SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo) |
172 : | maechler | 1598 | { |
173 : | maechler | 1960 | CHM_SP chx = AS_CHM_SP(x), chgx; |
174 : | maechler | 2113 | int uploT = (*CHAR(STRING_ELT(uplo,0)) == 'U') ? 1 : -1; |
175 : | maechler | 1736 | int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
176 : | maechler | 1960 | R_CheckStack(); |
177 : | maechler | 1598 | |
178 : | maechler | 1618 | chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c); |
179 : | maechler | 1598 | /* xtype: pattern, "real", complex or .. */ |
180 : | return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", | ||
181 : | GET_SLOT(x, Matrix_DimNamesSym)); | ||
182 : | } | ||
183 : | |||
184 : | bates | 1369 | SEXP Csparse_transpose(SEXP x, SEXP tri) |
185 : | bates | 922 | { |
186 : | maechler | 1921 | /* TODO: lgCMatrix & igC* currently go via double prec. cholmod - |
187 : | * since cholmod (& cs) lacks sparse 'int' matrices */ | ||
188 : | maechler | 1960 | CHM_SP chx = AS_CHM_SP(x); |
189 : | maechler | 1736 | int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
190 : | maechler | 1960 | CHM_SP chxt = cholmod_transpose(chx, chx->xtype, &c); |
191 : | bates | 1366 | SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp; |
192 : | bates | 1867 | int tr = asLogical(tri); |
193 : | maechler | 1960 | R_CheckStack(); |
194 : | bates | 1369 | |
195 : | bates | 1366 | tmp = VECTOR_ELT(dn, 0); /* swap the dimnames */ |
196 : | SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1)); | ||
197 : | SET_VECTOR_ELT(dn, 1, tmp); | ||
198 : | UNPROTECT(1); | ||
199 : | bates | 1867 | return chm_sparse_to_SEXP(chxt, 1, /* SWAP 'uplo' for triangular */ |
200 : | tr ? ((*uplo_P(x) == 'U') ? -1 : 1) : 0, | ||
201 : | Rkind, tr ? diag_P(x) : "", dn); | ||
202 : | bates | 922 | } |
203 : | |||
204 : | SEXP Csparse_Csparse_prod(SEXP a, SEXP b) | ||
205 : | { | ||
206 : | maechler | 2120 | CHM_SP |
207 : | cha = AS_CHM_SP(Csparse_diagU2N(a)), | ||
208 : | chb = AS_CHM_SP(Csparse_diagU2N(b)), | ||
209 : | chc = cholmod_ssmult(cha, chb, 0, cha->xtype, 1, &c); | ||
210 : | bates | 1366 | SEXP dn = allocVector(VECSXP, 2); |
211 : | maechler | 1960 | R_CheckStack(); |
212 : | bates | 922 | |
213 : | bates | 1366 | SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
214 : | duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); | ||
215 : | SET_VECTOR_ELT(dn, 1, | ||
216 : | duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1))); | ||
217 : | maechler | 1548 | return chm_sparse_to_SEXP(chc, 1, 0, 0, "", dn); |
218 : | bates | 922 | } |
219 : | |||
220 : | maechler | 1659 | SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b, SEXP trans) |
221 : | bates | 1657 | { |
222 : | maechler | 1659 | int tr = asLogical(trans); |
223 : | maechler | 2120 | CHM_SP |
224 : | cha = AS_CHM_SP(Csparse_diagU2N(a)), | ||
225 : | chb = AS_CHM_SP(Csparse_diagU2N(b)), | ||
226 : | chTr, chc; | ||
227 : | bates | 1657 | SEXP dn = allocVector(VECSXP, 2); |
228 : | maechler | 1960 | R_CheckStack(); |
229 : | bates | 1657 | |
230 : | maechler | 1960 | chTr = cholmod_transpose((tr) ? chb : cha, chb->xtype, &c); |
231 : | maechler | 1659 | chc = cholmod_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb, |
232 : | 0, cha->xtype, 1, &c); | ||
233 : | maechler | 1960 | cholmod_free_sparse(&chTr, &c); |
234 : | maechler | 1659 | |
235 : | bates | 1657 | SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
236 : | maechler | 1659 | duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), (tr) ? 0 : 1))); |
237 : | bates | 1657 | SET_VECTOR_ELT(dn, 1, |
238 : | maechler | 1659 | duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), (tr) ? 0 : 1))); |
239 : | bates | 1657 | return chm_sparse_to_SEXP(chc, 1, 0, 0, "", dn); |
240 : | } | ||
241 : | |||
242 : | bates | 922 | SEXP Csparse_dense_prod(SEXP a, SEXP b) |
243 : | { | ||
244 : | maechler | 2120 | CHM_SP cha = AS_CHM_SP(Csparse_diagU2N(a)); |
245 : | maechler | 1660 | SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
246 : | maechler | 1960 | CHM_DN chb = AS_CHM_DN(b_M); |
247 : | CHM_DN chc = cholmod_allocate_dense(cha->nrow, chb->ncol, cha->nrow, | ||
248 : | chb->xtype, &c); | ||
249 : | SEXP dn = PROTECT(allocVector(VECSXP, 2)); | ||
250 : | double one[] = {1,0}, zero[] = {0,0}; | ||
251 : | R_CheckStack(); | ||
252 : | bates | 922 | |
253 : | maechler | 1960 | cholmod_sdmult(cha, 0, one, zero, chb, chc, &c); |
254 : | maechler | 1660 | SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
255 : | duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); | ||
256 : | SET_VECTOR_ELT(dn, 1, | ||
257 : | duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1))); | ||
258 : | maechler | 1960 | UNPROTECT(2); |
259 : | maechler | 1660 | return chm_dense_to_SEXP(chc, 1, 0, dn); |
260 : | bates | 922 | } |
261 : | maechler | 925 | |
262 : | bates | 1067 | SEXP Csparse_dense_crossprod(SEXP a, SEXP b) |
263 : | { | ||
264 : | maechler | 2120 | CHM_SP cha = AS_CHM_SP(Csparse_diagU2N(a)); |
265 : | maechler | 1660 | SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
266 : | maechler | 1960 | CHM_DN chb = AS_CHM_DN(b_M); |
267 : | CHM_DN chc = cholmod_allocate_dense(cha->ncol, chb->ncol, cha->ncol, | ||
268 : | chb->xtype, &c); | ||
269 : | SEXP dn = PROTECT(allocVector(VECSXP, 2)); | ||
270 : | double one[] = {1,0}, zero[] = {0,0}; | ||
271 : | R_CheckStack(); | ||
272 : | bates | 1067 | |
273 : | maechler | 1960 | cholmod_sdmult(cha, 1, one, zero, chb, chc, &c); |
274 : | maechler | 1660 | SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
275 : | duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1))); | ||
276 : | SET_VECTOR_ELT(dn, 1, | ||
277 : | duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1))); | ||
278 : | maechler | 1960 | UNPROTECT(2); |
279 : | maechler | 1660 | return chm_dense_to_SEXP(chc, 1, 0, dn); |
280 : | bates | 1067 | } |
281 : | |||
282 : | maechler | 1659 | /* Computes x'x or x x' -- see Csparse_Csparse_crossprod above for x'y and x y' */ |
283 : | bates | 928 | SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet) |
284 : | bates | 922 | { |
285 : | maechler | 957 | int trip = asLogical(triplet), |
286 : | tr = asLogical(trans); /* gets reversed because _aat is tcrossprod */ | ||
287 : | maechler | 2120 | CHM_TR cht = trip ? AS_CHM_TR(Tsparse_diagU2N(x)) : (CHM_TR) NULL; |
288 : | maechler | 1960 | CHM_SP chcp, chxt, |
289 : | maechler | 2120 | chx = (trip ? |
290 : | cholmod_triplet_to_sparse(cht, cht->nnz, &c) : | ||
291 : | AS_CHM_SP(Csparse_diagU2N(x))); | ||
292 : | bates | 1366 | SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
293 : | maechler | 1960 | R_CheckStack(); |
294 : | bates | 922 | |
295 : | maechler | 1960 | if (!tr) chxt = cholmod_transpose(chx, chx->xtype, &c); |
296 : | bates | 928 | chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c); |
297 : | maechler | 2120 | if(!chcp) { |
298 : | UNPROTECT(1); | ||
299 : | error(_("Csparse_crossprod(): error return from cholmod_aat()")); | ||
300 : | } | ||
301 : | bates | 1360 | cholmod_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c); |
302 : | chcp->stype = 1; | ||
303 : | maechler | 1960 | if (trip) cholmod_free_sparse(&chx, &c); |
304 : | bates | 923 | if (!tr) cholmod_free_sparse(&chxt, &c); |
305 : | maechler | 1960 | SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
306 : | bates | 1366 | duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), |
307 : | maechler | 1660 | (tr) ? 0 : 1))); |
308 : | bates | 1366 | SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(dn, 0))); |
309 : | UNPROTECT(1); | ||
310 : | maechler | 1548 | return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn); |
311 : | bates | 922 | } |
312 : | bates | 923 | |
313 : | maechler | 1618 | SEXP Csparse_drop(SEXP x, SEXP tol) |
314 : | { | ||
315 : | maechler | 1960 | CHM_SP chx = AS_CHM_SP(x); |
316 : | CHM_SP ans = cholmod_copy(chx, chx->stype, chx->xtype, &c); | ||
317 : | maechler | 1618 | double dtol = asReal(tol); |
318 : | maechler | 1736 | int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
319 : | maechler | 1960 | R_CheckStack(); |
320 : | maechler | 1618 | |
321 : | if(!cholmod_drop(dtol, ans, &c)) | ||
322 : | error(_("cholmod_drop() failed")); | ||
323 : | maechler | 1736 | return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", |
324 : | GET_SLOT(x, Matrix_DimNamesSym)); | ||
325 : | maechler | 1618 | } |
326 : | |||
327 : | bates | 1218 | SEXP Csparse_horzcat(SEXP x, SEXP y) |
328 : | { | ||
329 : | maechler | 1960 | CHM_SP chx = AS_CHM_SP(x), chy = AS_CHM_SP(y); |
330 : | maechler | 1548 | int Rkind = 0; /* only for "d" - FIXME */ |
331 : | maechler | 1960 | R_CheckStack(); |
332 : | maechler | 1375 | |
333 : | bates | 1366 | /* FIXME: currently drops dimnames */ |
334 : | maechler | 1960 | return chm_sparse_to_SEXP(cholmod_horzcat(chx, chy, 1, &c), |
335 : | 1, 0, Rkind, "", R_NilValue); | ||
336 : | bates | 1218 | } |
337 : | |||
338 : | SEXP Csparse_vertcat(SEXP x, SEXP y) | ||
339 : | { | ||
340 : | maechler | 1960 | CHM_SP chx = AS_CHM_SP(x), chy = AS_CHM_SP(y); |
341 : | maechler | 1548 | int Rkind = 0; /* only for "d" - FIXME */ |
342 : | maechler | 1960 | R_CheckStack(); |
343 : | maechler | 1375 | |
344 : | bates | 1366 | /* FIXME: currently drops dimnames */ |
345 : | maechler | 1960 | return chm_sparse_to_SEXP(cholmod_vertcat(chx, chy, 1, &c), |
346 : | 1, 0, Rkind, "", R_NilValue); | ||
347 : | bates | 1218 | } |
348 : | bates | 1265 | |
349 : | SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2) | ||
350 : | { | ||
351 : | maechler | 1960 | CHM_SP chx = AS_CHM_SP(x); |
352 : | maechler | 1736 | int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
353 : | maechler | 1960 | CHM_SP ans = cholmod_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c); |
354 : | R_CheckStack(); | ||
355 : | bates | 1265 | |
356 : | maechler | 1736 | return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", |
357 : | GET_SLOT(x, Matrix_DimNamesSym)); | ||
358 : | bates | 1265 | } |
359 : | bates | 1366 | |
360 : | SEXP Csparse_diagU2N(SEXP x) | ||
361 : | { | ||
362 : | maechler | 2120 | const char *cl = class_P(x); |
363 : | /* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */ | ||
364 : | if (cl[1] != 't' || *diag_P(x) != 'U') { | ||
365 : | /* "trivially fast" when not triangular (<==> no 'diag' slot), or not *unit* triangular */ | ||
366 : | maechler | 1708 | return (x); |
367 : | } | ||
368 : | else { | ||
369 : | maechler | 1960 | CHM_SP chx = AS_CHM_SP(x); |
370 : | CHM_SP eye = cholmod_speye(chx->nrow, chx->ncol, chx->xtype, &c); | ||
371 : | maechler | 1708 | double one[] = {1, 0}; |
372 : | maechler | 1960 | CHM_SP ans = cholmod_add(chx, eye, one, one, TRUE, TRUE, &c); |
373 : | maechler | 1710 | int uploT = (*uplo_P(x) == 'U') ? 1 : -1; |
374 : | maechler | 1736 | int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
375 : | bates | 1366 | |
376 : | maechler | 1960 | R_CheckStack(); |
377 : | cholmod_free_sparse(&eye, &c); | ||
378 : | maechler | 1708 | return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N", |
379 : | maechler | 1736 | GET_SLOT(x, Matrix_DimNamesSym)); |
380 : | maechler | 1708 | } |
381 : | bates | 1366 | } |
382 : | |||
383 : | SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j) | ||
384 : | { | ||
385 : | maechler | 1960 | CHM_SP chx = AS_CHM_SP(x); |
386 : | bates | 1366 | int rsize = (isNull(i)) ? -1 : LENGTH(i), |
387 : | csize = (isNull(j)) ? -1 : LENGTH(j); | ||
388 : | maechler | 1736 | int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
389 : | maechler | 1960 | R_CheckStack(); |
390 : | bates | 1366 | |
391 : | if (rsize >= 0 && !isInteger(i)) | ||
392 : | error(_("Index i must be NULL or integer")); | ||
393 : | if (csize >= 0 && !isInteger(j)) | ||
394 : | error(_("Index j must be NULL or integer")); | ||
395 : | maechler | 1736 | |
396 : | bates | 1366 | return chm_sparse_to_SEXP(cholmod_submatrix(chx, INTEGER(i), rsize, |
397 : | maechler | 1375 | INTEGER(j), csize, |
398 : | bates | 1366 | TRUE, TRUE, &c), |
399 : | maechler | 1736 | 1, 0, Rkind, "", |
400 : | /* FIXME: drops dimnames */ R_NilValue); | ||
401 : | bates | 1366 | } |
402 : | bates | 2049 | |
403 : | SEXP Csparse_MatrixMarket(SEXP x, SEXP fname) | ||
404 : | { | ||
405 : | FILE *f = fopen(CHAR(asChar(fname)), "w"); | ||
406 : | |||
407 : | if (!f) | ||
408 : | error(_("failure to open file \"%s\" for writing"), | ||
409 : | CHAR(asChar(fname))); | ||
410 : | maechler | 2120 | if (!cholmod_write_sparse(f, AS_CHM_SP(Csparse_diagU2N(x)), |
411 : | (CHM_SP)NULL, (char*) NULL, &c)) | ||
412 : | bates | 2049 | error(_("cholmod_write_sparse returned error code")); |
413 : | fclose(f); | ||
414 : | return R_NilValue; | ||
415 : | } |
R-Forge@R-project.org | ViewVC Help |
Powered by ViewVC 1.0.0 |