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