7 |
/* NB: we do *NOT* check a potential 'x' slot here, at all */ |
/* 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, k, ncol, nrow, sorted, |
Rboolean sorted, strictly; |
11 |
|
int j, k, |
12 |
*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), |
*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), |
13 |
|
nrow = dims[0], |
14 |
|
ncol = dims[1], |
15 |
*xp = INTEGER(pslot), |
*xp = INTEGER(pslot), |
16 |
*xi = INTEGER(islot); |
*xi = INTEGER(islot); |
17 |
|
|
|
nrow = dims[0]; |
|
|
ncol = dims[1]; |
|
18 |
if (length(pslot) != dims[1] + 1) |
if (length(pslot) != dims[1] + 1) |
19 |
return mkString(_("slot p must have length = ncol(.) + 1")); |
return mkString(_("slot p must have length = ncol(.) + 1")); |
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 |
return |
24 |
mkString(_("last element of slot p must match length of slots i and x")); |
mkString(_("last element of slot p must match length of slots i and x")); |
25 |
for (j = 0; j < length(islot); j++) { |
for (j = 0; j < length(islot); j++) { |
26 |
if (xi[j] < 0 || xi[j] >= nrow) |
if (xi[j] < 0 || xi[j] >= nrow) |
27 |
return mkString(_("all row indices must be between 0 and nrow-1")); |
return mkString(_("all row indices must be between 0 and nrow-1")); |
28 |
} |
} |
29 |
sorted = TRUE; |
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 |
for (k = xp[j] + 1; k < xp[j + 1]; k++) |
if(sorted) |
34 |
if (xi[k] < xi[k - 1]) sorted = FALSE; |
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 |
} |
} |
41 |
if (!sorted) { |
if (!sorted) { |
42 |
cholmod_sparse *chx = as_cholmod_sparse(x); |
CHM_SP chx = AS_CHM_SP(x); |
43 |
|
R_CheckStack(); |
44 |
|
|
45 |
cholmod_sort(chx, &c); |
cholmod_sort(chx, &c); |
46 |
Free(chx); |
/* 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); |
58 |
|
} |
59 |
|
|
60 |
|
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); |
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) |
SEXP Csparse_to_dense(SEXP x) |
110 |
{ |
{ |
111 |
cholmod_sparse *chxs = as_cholmod_sparse(x); |
CHM_SP chxs = AS_CHM_SP(x); |
112 |
cholmod_dense *chxd = cholmod_sparse_to_dense(chxs, &c); |
/* 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 |
Free(chxs); |
return chm_dense_to_SEXP(chxd, 1, Rkind, GET_SLOT(x, Matrix_DimNamesSym)); |
|
return chm_dense_to_SEXP(chxd, 1, Real_kind(x)); |
|
120 |
} |
} |
121 |
|
|
122 |
SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri) |
SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri) |
123 |
{ |
{ |
124 |
cholmod_sparse *chxs = as_cholmod_sparse(x); |
CHM_SP chxs = AS_CHM_SP(x); |
125 |
cholmod_sparse |
CHM_SP chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c); |
126 |
*chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c); |
int tr = asLogical(tri); |
127 |
int uploT = 0; char *diag = ""; |
R_CheckStack(); |
128 |
|
|
129 |
Free(chxs); |
return chm_sparse_to_SEXP(chxcp, 1/*do_free*/, |
130 |
if (asLogical(tri)) { /* triangular sparse matrices */ |
tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, |
131 |
uploT = (strcmp(CHAR(asChar(GET_SLOT(x, Matrix_uploSym))), "U")) ? |
0, tr ? diag_P(x) : "", |
|
-1 : 1; |
|
|
diag = CHAR(asChar(GET_SLOT(x, Matrix_diagSym))); |
|
|
} |
|
|
return chm_sparse_to_SEXP(chxcp, 1, uploT, 0, diag, |
|
132 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
133 |
} |
} |
134 |
|
|
135 |
SEXP Csparse_to_matrix(SEXP x) |
SEXP Csparse_to_matrix(SEXP x) |
136 |
{ |
{ |
137 |
cholmod_sparse *chxs = as_cholmod_sparse(x); |
return chm_dense_to_matrix(cholmod_sparse_to_dense(AS_CHM_SP(x), &c), |
138 |
cholmod_dense *chxd = cholmod_sparse_to_dense(chxs, &c); |
1 /*do_free*/, GET_SLOT(x, Matrix_DimNamesSym)); |
|
|
|
|
Free(chxs); |
|
|
return chm_dense_to_matrix(chxd, 1, |
|
|
GET_SLOT(x, Matrix_DimNamesSym)); |
|
139 |
} |
} |
140 |
|
|
141 |
SEXP Csparse_to_Tsparse(SEXP x, SEXP tri) |
SEXP Csparse_to_Tsparse(SEXP x, SEXP tri) |
142 |
{ |
{ |
143 |
cholmod_sparse *chxs = as_cholmod_sparse(x); |
CHM_SP chxs = AS_CHM_SP(x); |
144 |
cholmod_triplet *chxt = cholmod_sparse_to_triplet(chxs, &c); |
CHM_TR chxt = cholmod_sparse_to_triplet(chxs, &c); |
145 |
int uploT = 0; |
int tr = asLogical(tri); |
146 |
char *diag = ""; |
int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
147 |
int Rkind = (chxs->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
R_CheckStack(); |
148 |
|
|
149 |
Free(chxs); |
return chm_triplet_to_SEXP(chxt, 1, |
150 |
if (asLogical(tri)) { /* triangular sparse matrices */ |
tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, |
151 |
uploT = (*uplo_P(x) == 'U') ? -1 : 1; |
Rkind, tr ? diag_P(x) : "", |
|
diag = diag_P(x); |
|
|
} |
|
|
return chm_triplet_to_SEXP(chxt, 1, uploT, Rkind, diag, |
|
152 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
153 |
} |
} |
154 |
|
|
155 |
/* this used to be called sCMatrix_to_gCMatrix(..) [in ./dsCMatrix.c ]: */ |
/* this used to be called sCMatrix_to_gCMatrix(..) [in ./dsCMatrix.c ]: */ |
156 |
SEXP Csparse_symmetric_to_general(SEXP x) |
SEXP Csparse_symmetric_to_general(SEXP x) |
157 |
{ |
{ |
158 |
cholmod_sparse *chx = as_cholmod_sparse(x), *chgx; |
CHM_SP chx = AS_CHM_SP(x), chgx; |
159 |
int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
160 |
|
R_CheckStack(); |
161 |
|
|
162 |
if (!(chx->stype)) |
if (!(chx->stype)) |
163 |
error(_("Nonsymmetric matrix in Csparse_symmetric_to_general")); |
error(_("Nonsymmetric matrix in Csparse_symmetric_to_general")); |
164 |
chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c); |
chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c); |
165 |
/* xtype: pattern, "real", complex or .. */ |
/* xtype: pattern, "real", complex or .. */ |
|
Free(chx); |
|
166 |
return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
167 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
168 |
} |
} |
169 |
|
|
170 |
SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo) |
SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo) |
171 |
{ |
{ |
172 |
cholmod_sparse *chx = as_cholmod_sparse(x), *chgx; |
CHM_SP chx = AS_CHM_SP(x), chgx; |
173 |
int uploT = (*CHAR(asChar(uplo)) == 'U') ? -1 : 1; |
int uploT = (*CHAR(asChar(uplo)) == 'U') ? 1 : -1; |
174 |
int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
175 |
|
R_CheckStack(); |
176 |
|
|
177 |
chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c); |
chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c); |
178 |
/* xtype: pattern, "real", complex or .. */ |
/* xtype: pattern, "real", complex or .. */ |
|
Free(chx); |
|
179 |
return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
180 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
181 |
} |
} |
182 |
|
|
183 |
SEXP Csparse_transpose(SEXP x, SEXP tri) |
SEXP Csparse_transpose(SEXP x, SEXP tri) |
184 |
{ |
{ |
185 |
cholmod_sparse *chx = as_cholmod_sparse(x); |
/* TODO: lgCMatrix & igC* currently go via double prec. cholmod - |
186 |
int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
* since cholmod (& cs) lacks sparse 'int' matrices */ |
187 |
cholmod_sparse *chxt = cholmod_transpose(chx, (int) chx->xtype, &c); |
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; |
SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp; |
191 |
int uploT = 0; char *diag = ""; |
int tr = asLogical(tri); |
192 |
|
R_CheckStack(); |
193 |
|
|
|
Free(chx); |
|
194 |
tmp = VECTOR_ELT(dn, 0); /* swap the dimnames */ |
tmp = VECTOR_ELT(dn, 0); /* swap the dimnames */ |
195 |
SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1)); |
SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1)); |
196 |
SET_VECTOR_ELT(dn, 1, tmp); |
SET_VECTOR_ELT(dn, 1, tmp); |
197 |
UNPROTECT(1); |
UNPROTECT(1); |
198 |
if (asLogical(tri)) { /* triangular sparse matrices */ |
return chm_sparse_to_SEXP(chxt, 1, /* SWAP 'uplo' for triangular */ |
199 |
uploT = (*uplo_P(x) == 'U') ? -1 : 1; |
tr ? ((*uplo_P(x) == 'U') ? -1 : 1) : 0, |
200 |
diag = diag_P(x); |
Rkind, tr ? diag_P(x) : "", dn); |
|
} |
|
|
return chm_sparse_to_SEXP(chxt, 1, uploT, Rkind, diag, dn); |
|
201 |
} |
} |
202 |
|
|
203 |
SEXP Csparse_Csparse_prod(SEXP a, SEXP b) |
SEXP Csparse_Csparse_prod(SEXP a, SEXP b) |
204 |
{ |
{ |
205 |
cholmod_sparse *cha = as_cholmod_sparse(a), |
CHM_SP cha = AS_CHM_SP(a), chb = AS_CHM_SP(b); |
206 |
*chb = as_cholmod_sparse(b); |
CHM_SP chc = cholmod_ssmult(cha, chb, 0, cha->xtype, 1, &c); |
|
cholmod_sparse *chc = cholmod_ssmult(cha, chb, 0, cha->xtype, 1, &c); |
|
207 |
SEXP dn = allocVector(VECSXP, 2); |
SEXP dn = allocVector(VECSXP, 2); |
208 |
|
R_CheckStack(); |
209 |
|
|
|
Free(cha); Free(chb); |
|
210 |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
211 |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); |
212 |
SET_VECTOR_ELT(dn, 1, |
SET_VECTOR_ELT(dn, 1, |
214 |
return chm_sparse_to_SEXP(chc, 1, 0, 0, "", dn); |
return chm_sparse_to_SEXP(chc, 1, 0, 0, "", dn); |
215 |
} |
} |
216 |
|
|
217 |
SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b) |
SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b, SEXP trans) |
218 |
{ |
{ |
219 |
cholmod_sparse *cha = as_cholmod_sparse(a), |
int tr = asLogical(trans); |
220 |
*chb = as_cholmod_sparse(b); |
CHM_SP cha = AS_CHM_SP(a), chb = AS_CHM_SP(b), chTr, chc; |
|
cholmod_sparse *chta = cholmod_transpose(cha, 1, &c); |
|
|
cholmod_sparse *chc = cholmod_ssmult(chta, chb, 0, cha->xtype, 1, &c); |
|
221 |
SEXP dn = allocVector(VECSXP, 2); |
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 |
|
|
|
Free(cha); Free(chb); cholmod_free_sparse(&chta, &c); |
|
229 |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
230 |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1))); |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), (tr) ? 0 : 1))); |
231 |
SET_VECTOR_ELT(dn, 1, |
SET_VECTOR_ELT(dn, 1, |
232 |
duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1))); |
duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), (tr) ? 0 : 1))); |
233 |
return chm_sparse_to_SEXP(chc, 1, 0, 0, "", dn); |
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 |
cholmod_sparse *cha = as_cholmod_sparse(a); |
CHM_SP cha = AS_CHM_SP(a); |
239 |
cholmod_dense *chb = as_cholmod_dense(PROTECT(mMatrix_as_dgeMatrix(b))); |
SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
240 |
cholmod_dense *chc = |
CHM_DN chb = AS_CHM_DN(b_M); |
241 |
cholmod_allocate_dense(cha->nrow, chb->ncol, cha->nrow, chb->xtype, &c); |
CHM_DN chc = cholmod_allocate_dense(cha->nrow, chb->ncol, cha->nrow, |
242 |
double alpha[] = {1,0}, beta[] = {0,0}; |
chb->xtype, &c); |
243 |
|
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
244 |
|
double one[] = {1,0}, zero[] = {0,0}; |
245 |
|
R_CheckStack(); |
246 |
|
|
247 |
cholmod_sdmult(cha, 0, alpha, beta, chb, chc, &c); |
cholmod_sdmult(cha, 0, one, zero, chb, chc, &c); |
248 |
Free(cha); Free(chb); |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
249 |
UNPROTECT(1); |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); |
250 |
return chm_dense_to_SEXP(chc, 1, 0); |
SET_VECTOR_ELT(dn, 1, |
251 |
|
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) |
SEXP Csparse_dense_crossprod(SEXP a, SEXP b) |
257 |
{ |
{ |
258 |
cholmod_sparse *cha = as_cholmod_sparse(a); |
CHM_SP cha = AS_CHM_SP(a); |
259 |
cholmod_dense *chb = as_cholmod_dense(PROTECT(mMatrix_as_dgeMatrix(b))); |
SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
260 |
cholmod_dense *chc = |
CHM_DN chb = AS_CHM_DN(b_M); |
261 |
cholmod_allocate_dense(cha->ncol, chb->ncol, cha->ncol, chb->xtype, &c); |
CHM_DN chc = cholmod_allocate_dense(cha->ncol, chb->ncol, cha->ncol, |
262 |
double alpha[] = {1,0}, beta[] = {0,0}; |
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, alpha, beta, chb, chc, &c); |
cholmod_sdmult(cha, 1, one, zero, chb, chc, &c); |
268 |
Free(cha); Free(chb); |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
269 |
UNPROTECT(1); |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1))); |
270 |
return chm_dense_to_SEXP(chc, 1, 0); |
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 |
{ |
{ |
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); |
|
*chx = trip ? cholmod_triplet_to_sparse(cht, cht->nnz, &c) |
|
|
: as_cholmod_sparse(x); |
|
284 |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
285 |
|
R_CheckStack(); |
286 |
|
|
287 |
if (!tr) |
if (!tr) chxt = cholmod_transpose(chx, chx->xtype, &c); |
|
chxt = cholmod_transpose(chx, 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()")); |
|
error(_("Csparse_crossprod(): error return from cholmod_aat()")); |
|
290 |
cholmod_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c); |
cholmod_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c); |
291 |
chcp->stype = 1; |
chcp->stype = 1; |
292 |
if (trip) { |
if (trip) cholmod_free_sparse(&chx, &c); |
|
cholmod_free_sparse(&chx, &c); |
|
|
Free(cht); |
|
|
} else { |
|
|
Free(chx); |
|
|
} |
|
293 |
if (!tr) cholmod_free_sparse(&chxt, &c); |
if (!tr) cholmod_free_sparse(&chxt, &c); |
294 |
/* create dimnames */ |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
|
SET_VECTOR_ELT(dn, 0, |
|
295 |
duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), |
duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), |
296 |
(tr) ? 1 : 0))); |
(tr) ? 0 : 1))); |
297 |
SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(dn, 0))); |
SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(dn, 0))); |
298 |
UNPROTECT(1); |
UNPROTECT(1); |
299 |
return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn); |
return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn); |
301 |
|
|
302 |
SEXP Csparse_drop(SEXP x, SEXP tol) |
SEXP Csparse_drop(SEXP x, SEXP tol) |
303 |
{ |
{ |
304 |
cholmod_sparse *chx = as_cholmod_sparse(x), |
CHM_SP chx = AS_CHM_SP(x); |
305 |
*ans = cholmod_copy(chx, chx->stype, chx->xtype, &c); |
CHM_SP ans = cholmod_copy(chx, chx->stype, chx->xtype, &c); |
306 |
double dtol = asReal(tol); |
double dtol = asReal(tol); |
307 |
int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
308 |
|
R_CheckStack(); |
309 |
|
|
310 |
if(!cholmod_drop(dtol, ans, &c)) |
if(!cholmod_drop(dtol, ans, &c)) |
311 |
error(_("cholmod_drop() failed")); |
error(_("cholmod_drop() failed")); |
312 |
Free(chx); |
return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", |
313 |
/* FIXME: currently drops dimnames */ |
GET_SLOT(x, Matrix_DimNamesSym)); |
|
return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue); |
|
314 |
} |
} |
315 |
|
|
|
|
|
316 |
SEXP Csparse_horzcat(SEXP x, SEXP y) |
SEXP Csparse_horzcat(SEXP x, SEXP y) |
317 |
{ |
{ |
318 |
cholmod_sparse *chx = as_cholmod_sparse(x), |
CHM_SP chx = AS_CHM_SP(x), chy = AS_CHM_SP(y); |
|
*chy = as_cholmod_sparse(y), *ans; |
|
319 |
int Rkind = 0; /* only for "d" - FIXME */ |
int Rkind = 0; /* only for "d" - FIXME */ |
320 |
|
R_CheckStack(); |
321 |
|
|
|
ans = cholmod_horzcat(chx, chy, 1, &c); |
|
|
Free(chx); Free(chy); |
|
322 |
/* FIXME: currently drops dimnames */ |
/* FIXME: currently drops dimnames */ |
323 |
return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue); |
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) |
SEXP Csparse_vertcat(SEXP x, SEXP y) |
328 |
{ |
{ |
329 |
cholmod_sparse *chx = as_cholmod_sparse(x), |
CHM_SP chx = AS_CHM_SP(x), chy = AS_CHM_SP(y); |
|
*chy = as_cholmod_sparse(y), *ans; |
|
330 |
int Rkind = 0; /* only for "d" - FIXME */ |
int Rkind = 0; /* only for "d" - FIXME */ |
331 |
|
R_CheckStack(); |
332 |
|
|
|
ans = cholmod_vertcat(chx, chy, 1, &c); |
|
|
Free(chx); Free(chy); |
|
333 |
/* FIXME: currently drops dimnames */ |
/* FIXME: currently drops dimnames */ |
334 |
return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue); |
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) |
SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2) |
339 |
{ |
{ |
340 |
cholmod_sparse *chx = as_cholmod_sparse(x), *ans; |
CHM_SP chx = AS_CHM_SP(x); |
341 |
int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
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 |
ans = cholmod_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c); |
return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", |
346 |
Free(chx); |
GET_SLOT(x, Matrix_DimNamesSym)); |
|
return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue); |
|
347 |
} |
} |
348 |
|
|
349 |
SEXP Csparse_diagU2N(SEXP x) |
SEXP Csparse_diagU2N(SEXP x) |
350 |
{ |
{ |
351 |
cholmod_sparse *chx = as_cholmod_sparse(x); |
if (*diag_P(x) != 'U') {/* "trivially fast" when there's no 'diag' slot at all */ |
352 |
cholmod_sparse *eye = cholmod_speye(chx->nrow, chx->ncol, chx->xtype, &c); |
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}; |
double one[] = {1, 0}; |
358 |
cholmod_sparse *ans = cholmod_add(chx, eye, one, one, TRUE, TRUE, &c); |
CHM_SP ans = cholmod_add(chx, eye, one, one, TRUE, TRUE, &c); |
359 |
int uploT = (strcmp(CHAR(asChar(GET_SLOT(x, Matrix_uploSym))), "U")) ? |
int uploT = (*uplo_P(x) == 'U') ? 1 : -1; |
360 |
-1 : 1; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
|
int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
|
361 |
|
|
362 |
Free(chx); cholmod_free_sparse(&eye, &c); |
R_CheckStack(); |
363 |
|
cholmod_free_sparse(&eye, &c); |
364 |
return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N", |
return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N", |
365 |
duplicate(GET_SLOT(x, Matrix_DimNamesSym))); |
GET_SLOT(x, Matrix_DimNamesSym)); |
366 |
|
} |
367 |
} |
} |
368 |
|
|
369 |
SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j) |
SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j) |
370 |
{ |
{ |
371 |
cholmod_sparse *chx = as_cholmod_sparse(x); |
CHM_SP chx = AS_CHM_SP(x); |
372 |
int rsize = (isNull(i)) ? -1 : LENGTH(i), |
int rsize = (isNull(i)) ? -1 : LENGTH(i), |
373 |
csize = (isNull(j)) ? -1 : LENGTH(j); |
csize = (isNull(j)) ? -1 : LENGTH(j); |
374 |
int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
375 |
|
R_CheckStack(); |
376 |
|
|
377 |
if (rsize >= 0 && !isInteger(i)) |
if (rsize >= 0 && !isInteger(i)) |
378 |
error(_("Index i must be NULL or integer")); |
error(_("Index i must be NULL or integer")); |
379 |
if (csize >= 0 && !isInteger(j)) |
if (csize >= 0 && !isInteger(j)) |
380 |
error(_("Index j must be NULL or integer")); |
error(_("Index j must be NULL or integer")); |
381 |
|
|
382 |
return chm_sparse_to_SEXP(cholmod_submatrix(chx, INTEGER(i), rsize, |
return chm_sparse_to_SEXP(cholmod_submatrix(chx, INTEGER(i), rsize, |
383 |
INTEGER(j), csize, |
INTEGER(j), csize, |
384 |
TRUE, TRUE, &c), |
TRUE, TRUE, &c), |
385 |
1, 0, Rkind, "", R_NilValue); |
1, 0, Rkind, "", |
386 |
|
/* FIXME: drops dimnames */ R_NilValue); |
387 |
|
} |
388 |
|
|
389 |
|
SEXP Csparse_MatrixMarket(SEXP x, SEXP fname) |
390 |
|
{ |
391 |
|
FILE *f = fopen(CHAR(asChar(fname)), "w"); |
392 |
|
|
393 |
|
if (!f) |
394 |
|
error(_("failure to open file \"%s\" for writing"), |
395 |
|
CHAR(asChar(fname))); |
396 |
|
if (!cholmod_write_sparse(f, AS_CHM_SP(x), (CHM_SP)NULL, |
397 |
|
(char*) NULL, &c)) |
398 |
|
error(_("cholmod_write_sparse returned error code")); |
399 |
|
fclose(f); |
400 |
|
return R_NilValue; |
401 |
} |
} |