1 |
/* Sparse matrices in compressed column-oriented form */ |
/* Sparse matrices in compressed column-oriented form */ |
2 |
#include "Csparse.h" |
#include "Csparse.h" |
3 |
|
#include "Tsparse.h" |
4 |
#include "chm_common.h" |
#include "chm_common.h" |
5 |
|
|
6 |
SEXP Csparse_validate(SEXP x) |
/** "Cheap" C version of Csparse_validate() - *not* sorting : */ |
7 |
|
Rboolean isValid_Csparse(SEXP x) |
8 |
{ |
{ |
9 |
/* NB: we do *NOT* check a potential 'x' slot here, at all */ |
/* NB: we do *NOT* check a potential 'x' slot here, at all */ |
10 |
SEXP pslot = GET_SLOT(x, Matrix_pSym), |
SEXP pslot = GET_SLOT(x, Matrix_pSym), |
11 |
islot = GET_SLOT(x, Matrix_iSym); |
islot = GET_SLOT(x, Matrix_iSym); |
12 |
int j, k, ncol, nrow, sorted, |
int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), j, |
13 |
|
nrow = dims[0], |
14 |
|
ncol = dims[1], |
15 |
|
*xp = INTEGER(pslot), |
16 |
|
*xi = INTEGER(islot); |
17 |
|
|
18 |
|
if (length(pslot) != dims[1] + 1) |
19 |
|
return FALSE; |
20 |
|
if (xp[0] != 0) |
21 |
|
return FALSE; |
22 |
|
if (length(islot) < xp[ncol]) /* allow larger slots from over-allocation!*/ |
23 |
|
return FALSE; |
24 |
|
for (j = 0; j < xp[ncol]; j++) { |
25 |
|
if (xi[j] < 0 || xi[j] >= nrow) |
26 |
|
return FALSE; |
27 |
|
} |
28 |
|
for (j = 0; j < ncol; j++) { |
29 |
|
if (xp[j] > xp[j + 1]) |
30 |
|
return FALSE; |
31 |
|
} |
32 |
|
return TRUE; |
33 |
|
} |
34 |
|
|
35 |
|
SEXP Csparse_validate(SEXP x) { |
36 |
|
return Csparse_validate_(x, FALSE); |
37 |
|
} |
38 |
|
|
39 |
|
SEXP Csparse_validate2(SEXP x, SEXP maybe_modify) { |
40 |
|
return Csparse_validate_(x, asLogical(maybe_modify)); |
41 |
|
} |
42 |
|
|
43 |
|
SEXP Csparse_validate_(SEXP x, Rboolean maybe_modify) |
44 |
|
{ |
45 |
|
/* NB: we do *NOT* check a potential 'x' slot here, at all */ |
46 |
|
SEXP pslot = GET_SLOT(x, Matrix_pSym), |
47 |
|
islot = GET_SLOT(x, Matrix_iSym); |
48 |
|
Rboolean sorted, strictly; |
49 |
|
int j, k, |
50 |
*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), |
*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), |
51 |
|
nrow = dims[0], |
52 |
|
ncol = dims[1], |
53 |
*xp = INTEGER(pslot), |
*xp = INTEGER(pslot), |
54 |
*xi = INTEGER(islot); |
*xi = INTEGER(islot); |
55 |
|
|
|
nrow = dims[0]; |
|
|
ncol = dims[1]; |
|
56 |
if (length(pslot) != dims[1] + 1) |
if (length(pslot) != dims[1] + 1) |
57 |
return mkString(_("slot p must have length = ncol(.) + 1")); |
return mkString(_("slot p must have length = ncol(.) + 1")); |
58 |
if (xp[0] != 0) |
if (xp[0] != 0) |
59 |
return mkString(_("first element of slot p must be zero")); |
return mkString(_("first element of slot p must be zero")); |
60 |
if (length(islot) != xp[ncol]) |
if (length(islot) < xp[ncol]) /* allow larger slots from over-allocation!*/ |
61 |
return |
return |
62 |
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")); |
63 |
for (j = 0; j < length(islot); j++) { |
for (j = 0; j < xp[ncol]; j++) { |
64 |
if (xi[j] < 0 || xi[j] >= nrow) |
if (xi[j] < 0 || xi[j] >= nrow) |
65 |
return mkString(_("all row indices must be between 0 and nrow-1")); |
return mkString(_("all row indices must be between 0 and nrow-1")); |
66 |
} |
} |
67 |
sorted = TRUE; |
sorted = TRUE; strictly = TRUE; |
68 |
for (j = 0; j < ncol; j++) { |
for (j = 0; j < ncol; j++) { |
69 |
if (xp[j] > xp[j+1]) |
if (xp[j] > xp[j+1]) |
70 |
return mkString(_("slot p must be non-decreasing")); |
return mkString(_("slot p must be non-decreasing")); |
71 |
for (k = xp[j] + 1; k < xp[j + 1]; k++) |
if(sorted) /* only act if >= 2 entries in column j : */ |
72 |
if (xi[k] < xi[k - 1]) sorted = FALSE; |
for (k = xp[j] + 1; k < xp[j + 1]; k++) { |
73 |
|
if (xi[k] < xi[k - 1]) |
74 |
|
sorted = FALSE; |
75 |
|
else if (xi[k] == xi[k - 1]) |
76 |
|
strictly = FALSE; |
77 |
|
} |
78 |
} |
} |
79 |
if (!sorted) { |
if (!sorted) { |
80 |
cholmod_sparse *chx = as_cholmod_sparse(x); |
if(maybe_modify) { |
81 |
cholmod_sort(chx, &c); |
CHM_SP chx = (CHM_SP) alloca(sizeof(cholmod_sparse)); |
82 |
Free(chx); |
R_CheckStack(); |
83 |
|
as_cholmod_sparse(chx, x, FALSE, TRUE);/*-> cholmod_l_sort() ! */ |
84 |
|
/* as chx = AS_CHM_SP__(x) but ^^^^ sorting x in_place !!! */ |
85 |
|
|
86 |
|
/* Now re-check that row indices are *strictly* increasing |
87 |
|
* (and not just increasing) within each column : */ |
88 |
|
for (j = 0; j < ncol; j++) { |
89 |
|
for (k = xp[j] + 1; k < xp[j + 1]; k++) |
90 |
|
if (xi[k] == xi[k - 1]) |
91 |
|
return mkString(_("slot i is not *strictly* increasing inside a column (even after cholmod_l_sort)")); |
92 |
|
} |
93 |
|
} else { /* no modifying sorting : */ |
94 |
|
return mkString(_("row indices are not sorted within columns")); |
95 |
|
} |
96 |
|
} else if(!strictly) { /* sorted, but not strictly */ |
97 |
|
return mkString(_("slot i is not *strictly* increasing inside a column")); |
98 |
} |
} |
99 |
return ScalarLogical(1); |
return ScalarLogical(1); |
100 |
} |
} |
101 |
|
|
102 |
|
SEXP Rsparse_validate(SEXP x) |
103 |
|
{ |
104 |
|
/* NB: we do *NOT* check a potential 'x' slot here, at all */ |
105 |
|
SEXP pslot = GET_SLOT(x, Matrix_pSym), |
106 |
|
jslot = GET_SLOT(x, Matrix_jSym); |
107 |
|
Rboolean sorted, strictly; |
108 |
|
int i, k, |
109 |
|
*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), |
110 |
|
nrow = dims[0], |
111 |
|
ncol = dims[1], |
112 |
|
*xp = INTEGER(pslot), |
113 |
|
*xj = INTEGER(jslot); |
114 |
|
|
115 |
|
if (length(pslot) != dims[0] + 1) |
116 |
|
return mkString(_("slot p must have length = nrow(.) + 1")); |
117 |
|
if (xp[0] != 0) |
118 |
|
return mkString(_("first element of slot p must be zero")); |
119 |
|
if (length(jslot) < xp[nrow]) /* allow larger slots from over-allocation!*/ |
120 |
|
return |
121 |
|
mkString(_("last element of slot p must match length of slots j and x")); |
122 |
|
for (i = 0; i < length(jslot); i++) { |
123 |
|
if (xj[i] < 0 || xj[i] >= ncol) |
124 |
|
return mkString(_("all column indices must be between 0 and ncol-1")); |
125 |
|
} |
126 |
|
sorted = TRUE; strictly = TRUE; |
127 |
|
for (i = 0; i < nrow; i++) { |
128 |
|
if (xp[i] > xp[i+1]) |
129 |
|
return mkString(_("slot p must be non-decreasing")); |
130 |
|
if(sorted) |
131 |
|
for (k = xp[i] + 1; k < xp[i + 1]; k++) { |
132 |
|
if (xj[k] < xj[k - 1]) |
133 |
|
sorted = FALSE; |
134 |
|
else if (xj[k] == xj[k - 1]) |
135 |
|
strictly = FALSE; |
136 |
|
} |
137 |
|
} |
138 |
|
if (!sorted) |
139 |
|
/* cannot easily use cholmod_l_sort(.) ... -> "error out" :*/ |
140 |
|
return mkString(_("slot j is not increasing inside a column")); |
141 |
|
else if(!strictly) /* sorted, but not strictly */ |
142 |
|
return mkString(_("slot j is not *strictly* increasing inside a column")); |
143 |
|
|
144 |
|
return ScalarLogical(1); |
145 |
|
} |
146 |
|
|
147 |
|
|
148 |
|
/* Called from ../R/Csparse.R : */ |
149 |
|
/* Can only return [dln]geMatrix (no symm/triang); |
150 |
|
* FIXME: replace by non-CHOLMOD code ! */ |
151 |
SEXP Csparse_to_dense(SEXP x) |
SEXP Csparse_to_dense(SEXP x) |
152 |
{ |
{ |
153 |
cholmod_sparse *chxs = as_cholmod_sparse(x); |
CHM_SP chxs = AS_CHM_SP__(x); |
154 |
cholmod_dense *chxd = cholmod_sparse_to_dense(chxs, &c); |
/* This loses the symmetry property, since cholmod_dense has none, |
155 |
|
* BUT, much worse (FIXME!), it also transforms CHOLMOD_PATTERN ("n") matrices |
156 |
|
* to numeric (CHOLMOD_REAL) ones : */ |
157 |
|
CHM_DN chxd = cholmod_l_sparse_to_dense(chxs, &c); |
158 |
|
int Rkind = (chxs->xtype == CHOLMOD_PATTERN)? -1 : Real_kind(x); |
159 |
|
R_CheckStack(); |
160 |
|
|
161 |
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)); |
|
162 |
} |
} |
163 |
|
|
164 |
SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri) |
SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri) |
165 |
{ |
{ |
166 |
cholmod_sparse *chxs = as_cholmod_sparse(x); |
CHM_SP chxs = AS_CHM_SP__(x); |
167 |
cholmod_sparse |
CHM_SP chxcp = cholmod_l_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c); |
168 |
*chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c); |
int tr = asLogical(tri); |
169 |
int uploT = 0; char *diag = ""; |
R_CheckStack(); |
170 |
|
|
171 |
Free(chxs); |
return chm_sparse_to_SEXP(chxcp, 1/*do_free*/, |
172 |
if (asLogical(tri)) { /* triangular sparse matrices */ |
tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, |
173 |
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, |
|
174 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
175 |
} |
} |
176 |
|
|
177 |
SEXP Csparse_to_matrix(SEXP x) |
SEXP Csparse_to_matrix(SEXP x) |
178 |
{ |
{ |
179 |
cholmod_sparse *chxs = as_cholmod_sparse(x); |
return chm_dense_to_matrix(cholmod_l_sparse_to_dense(AS_CHM_SP__(x), &c), |
180 |
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)); |
|
181 |
} |
} |
182 |
|
|
183 |
SEXP Csparse_to_Tsparse(SEXP x, SEXP tri) |
SEXP Csparse_to_Tsparse(SEXP x, SEXP tri) |
184 |
{ |
{ |
185 |
cholmod_sparse *chxs = as_cholmod_sparse(x); |
CHM_SP chxs = AS_CHM_SP__(x); |
186 |
cholmod_triplet *chxt = cholmod_sparse_to_triplet(chxs, &c); |
CHM_TR chxt = cholmod_l_sparse_to_triplet(chxs, &c); |
187 |
int uploT = 0; |
int tr = asLogical(tri); |
188 |
char *diag = ""; |
int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
189 |
int Rkind = (chxs->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
R_CheckStack(); |
190 |
|
|
191 |
Free(chxs); |
return chm_triplet_to_SEXP(chxt, 1, |
192 |
if (asLogical(tri)) { /* triangular sparse matrices */ |
tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, |
193 |
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, |
|
194 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
195 |
} |
} |
196 |
|
|
197 |
/* this used to be called sCMatrix_to_gCMatrix(..) [in ./dsCMatrix.c ]: */ |
/* this used to be called sCMatrix_to_gCMatrix(..) [in ./dsCMatrix.c ]: */ |
198 |
SEXP Csparse_symmetric_to_general(SEXP x) |
SEXP Csparse_symmetric_to_general(SEXP x) |
199 |
{ |
{ |
200 |
cholmod_sparse *chx = as_cholmod_sparse(x), *chgx; |
CHM_SP chx = AS_CHM_SP__(x), chgx; |
201 |
int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
202 |
|
R_CheckStack(); |
203 |
|
|
204 |
if (!(chx->stype)) |
if (!(chx->stype)) |
205 |
error(_("Nonsymmetric matrix in Csparse_symmetric_to_general")); |
error(_("Nonsymmetric matrix in Csparse_symmetric_to_general")); |
206 |
chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c); |
chgx = cholmod_l_copy(chx, /* stype: */ 0, chx->xtype, &c); |
207 |
/* xtype: pattern, "real", complex or .. */ |
/* xtype: pattern, "real", complex or .. */ |
|
Free(chx); |
|
208 |
return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
209 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
210 |
} |
} |
211 |
|
|
212 |
SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo) |
SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo) |
213 |
{ |
{ |
214 |
cholmod_sparse *chx = as_cholmod_sparse(x), *chgx; |
CHM_SP chx = AS_CHM_SP__(x), chgx; |
215 |
int uploT = (*CHAR(asChar(uplo)) == 'U') ? -1 : 1; |
int uploT = (*CHAR(STRING_ELT(uplo,0)) == 'U') ? 1 : -1; |
216 |
int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
217 |
|
R_CheckStack(); |
218 |
|
|
219 |
chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c); |
chgx = cholmod_l_copy(chx, /* stype: */ uploT, chx->xtype, &c); |
220 |
/* xtype: pattern, "real", complex or .. */ |
/* xtype: pattern, "real", complex or .. */ |
|
Free(chx); |
|
221 |
return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
222 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
223 |
} |
} |
224 |
|
|
225 |
SEXP Csparse_transpose(SEXP x, SEXP tri) |
SEXP Csparse_transpose(SEXP x, SEXP tri) |
226 |
{ |
{ |
227 |
cholmod_sparse *chx = as_cholmod_sparse(x); |
/* TODO: lgCMatrix & igC* currently go via double prec. cholmod - |
228 |
int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
* since cholmod (& cs) lacks sparse 'int' matrices */ |
229 |
cholmod_sparse *chxt = cholmod_transpose(chx, (int) chx->xtype, &c); |
CHM_SP chx = AS_CHM_SP__(x); |
230 |
|
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
231 |
|
CHM_SP chxt = cholmod_l_transpose(chx, chx->xtype, &c); |
232 |
SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp; |
SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp; |
233 |
int uploT = 0; char *diag = ""; |
int tr = asLogical(tri); |
234 |
|
R_CheckStack(); |
235 |
|
|
|
Free(chx); |
|
236 |
tmp = VECTOR_ELT(dn, 0); /* swap the dimnames */ |
tmp = VECTOR_ELT(dn, 0); /* swap the dimnames */ |
237 |
SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1)); |
SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1)); |
238 |
SET_VECTOR_ELT(dn, 1, tmp); |
SET_VECTOR_ELT(dn, 1, tmp); |
239 |
UNPROTECT(1); |
UNPROTECT(1); |
240 |
if (asLogical(tri)) { /* triangular sparse matrices */ |
return chm_sparse_to_SEXP(chxt, 1, /* SWAP 'uplo' for triangular */ |
241 |
uploT = (*uplo_P(x) == 'U') ? -1 : 1; |
tr ? ((*uplo_P(x) == 'U') ? -1 : 1) : 0, |
242 |
diag = diag_P(x); |
Rkind, tr ? diag_P(x) : "", dn); |
|
} |
|
|
return chm_sparse_to_SEXP(chxt, 1, uploT, Rkind, diag, dn); |
|
243 |
} |
} |
244 |
|
|
245 |
SEXP Csparse_Csparse_prod(SEXP a, SEXP b) |
SEXP Csparse_Csparse_prod(SEXP a, SEXP b) |
246 |
{ |
{ |
247 |
cholmod_sparse *cha = as_cholmod_sparse(a), |
CHM_SP |
248 |
*chb = as_cholmod_sparse(b); |
cha = AS_CHM_SP(a), |
249 |
cholmod_sparse *chc = cholmod_ssmult(cha, chb, 0, cha->xtype, 1, &c); |
chb = AS_CHM_SP(b), |
250 |
|
chc = cholmod_l_ssmult(cha, chb, /*out_stype:*/ 0, |
251 |
|
/* values:= is_numeric (T/F) */ cha->xtype > 0, |
252 |
|
/*out sorted:*/ 1, &c); |
253 |
|
const char *cl_a = class_P(a), *cl_b = class_P(b); |
254 |
|
char diag[] = {'\0', '\0'}; |
255 |
|
int uploT = 0; |
256 |
SEXP dn = allocVector(VECSXP, 2); |
SEXP dn = allocVector(VECSXP, 2); |
257 |
|
R_CheckStack(); |
258 |
|
|
259 |
Free(cha); Free(chb); |
#ifdef DEBUG_Matrix |
260 |
|
Rprintf("DBG Csparse_C*_prod(%s, %s)\n", cl_a, cl_b); |
261 |
|
#endif |
262 |
|
|
263 |
|
/* Preserve triangularity and even unit-triangularity if appropriate. |
264 |
|
* Note that in that case, the multiplication itself should happen |
265 |
|
* faster. But there's no support for that in CHOLMOD */ |
266 |
|
|
267 |
|
/* UGLY hack -- rather should have (fast!) C-level version of |
268 |
|
* is(a, "triangularMatrix") etc */ |
269 |
|
if (cl_a[1] == 't' && cl_b[1] == 't') |
270 |
|
/* FIXME: fails for "Cholesky","BunchKaufmann"..*/ |
271 |
|
if(*uplo_P(a) == *uplo_P(b)) { /* both upper, or both lower tri. */ |
272 |
|
uploT = (*uplo_P(a) == 'U') ? 1 : -1; |
273 |
|
if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */ |
274 |
|
/* "remove the diagonal entries": */ |
275 |
|
chm_diagN2U(chc, uploT, /* do_realloc */ FALSE); |
276 |
|
diag[0]= 'U'; |
277 |
|
} |
278 |
|
else diag[0]= 'N'; |
279 |
|
} |
280 |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
281 |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); |
282 |
SET_VECTOR_ELT(dn, 1, |
SET_VECTOR_ELT(dn, 1, |
283 |
duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1))); |
duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1))); |
284 |
return chm_sparse_to_SEXP(chc, 1, 0, 0, "", dn); |
return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn); |
285 |
} |
} |
286 |
|
|
287 |
SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b) |
SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b, SEXP trans) |
288 |
{ |
{ |
289 |
cholmod_sparse *cha = as_cholmod_sparse(a), |
int tr = asLogical(trans); |
290 |
*chb = as_cholmod_sparse(b); |
CHM_SP |
291 |
cholmod_sparse *chta = cholmod_transpose(cha, 1, &c); |
cha = AS_CHM_SP(a), |
292 |
cholmod_sparse *chc = cholmod_ssmult(chta, chb, 0, cha->xtype, 1, &c); |
chb = AS_CHM_SP(b), |
293 |
|
chTr, chc; |
294 |
|
const char *cl_a = class_P(a), *cl_b = class_P(b); |
295 |
|
char diag[] = {'\0', '\0'}; |
296 |
|
int uploT = 0; |
297 |
SEXP dn = allocVector(VECSXP, 2); |
SEXP dn = allocVector(VECSXP, 2); |
298 |
|
R_CheckStack(); |
299 |
|
|
300 |
|
chTr = cholmod_l_transpose((tr) ? chb : cha, chb->xtype, &c); |
301 |
|
chc = cholmod_l_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb, |
302 |
|
/*out_stype:*/ 0, cha->xtype, /*out sorted:*/ 1, &c); |
303 |
|
cholmod_l_free_sparse(&chTr, &c); |
304 |
|
|
305 |
|
/* Preserve triangularity and unit-triangularity if appropriate; |
306 |
|
* see Csparse_Csparse_prod() for comments */ |
307 |
|
if (cl_a[1] == 't' && cl_b[1] == 't') |
308 |
|
if(*uplo_P(a) != *uplo_P(b)) { /* one 'U', the other 'L' */ |
309 |
|
uploT = (*uplo_P(b) == 'U') ? 1 : -1; |
310 |
|
if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */ |
311 |
|
chm_diagN2U(chc, uploT, /* do_realloc */ FALSE); |
312 |
|
diag[0]= 'U'; |
313 |
|
} |
314 |
|
else diag[0]= 'N'; |
315 |
|
} |
316 |
|
|
|
Free(cha); Free(chb); cholmod_free_sparse(&chta, &c); |
|
317 |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
318 |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1))); |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), (tr) ? 0 : 1))); |
319 |
SET_VECTOR_ELT(dn, 1, |
SET_VECTOR_ELT(dn, 1, |
320 |
duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1))); |
duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), (tr) ? 0 : 1))); |
321 |
return chm_sparse_to_SEXP(chc, 1, 0, 0, "", dn); |
return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn); |
322 |
} |
} |
323 |
|
|
324 |
SEXP Csparse_dense_prod(SEXP a, SEXP b) |
SEXP Csparse_dense_prod(SEXP a, SEXP b) |
325 |
{ |
{ |
326 |
cholmod_sparse *cha = as_cholmod_sparse(a); |
CHM_SP cha = AS_CHM_SP(a); |
327 |
cholmod_dense *chb = as_cholmod_dense(PROTECT(mMatrix_as_dgeMatrix(b))); |
SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
328 |
cholmod_dense *chc = |
CHM_DN chb = AS_CHM_DN(b_M); |
329 |
cholmod_allocate_dense(cha->nrow, chb->ncol, cha->nrow, chb->xtype, &c); |
CHM_DN chc = cholmod_l_allocate_dense(cha->nrow, chb->ncol, cha->nrow, |
330 |
double alpha[] = {1,0}, beta[] = {0,0}; |
chb->xtype, &c); |
331 |
|
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
332 |
|
double one[] = {1,0}, zero[] = {0,0}; |
333 |
|
R_CheckStack(); |
334 |
|
|
335 |
cholmod_sdmult(cha, 0, alpha, beta, chb, chc, &c); |
cholmod_l_sdmult(cha, 0, one, zero, chb, chc, &c); |
336 |
Free(cha); Free(chb); |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
337 |
UNPROTECT(1); |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); |
338 |
return chm_dense_to_SEXP(chc, 1, 0); |
SET_VECTOR_ELT(dn, 1, |
339 |
|
duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1))); |
340 |
|
UNPROTECT(2); |
341 |
|
return chm_dense_to_SEXP(chc, 1, 0, dn); |
342 |
} |
} |
343 |
|
|
344 |
SEXP Csparse_dense_crossprod(SEXP a, SEXP b) |
SEXP Csparse_dense_crossprod(SEXP a, SEXP b) |
345 |
{ |
{ |
346 |
cholmod_sparse *cha = as_cholmod_sparse(a); |
CHM_SP cha = AS_CHM_SP(a); |
347 |
cholmod_dense *chb = as_cholmod_dense(PROTECT(mMatrix_as_dgeMatrix(b))); |
SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
348 |
cholmod_dense *chc = |
CHM_DN chb = AS_CHM_DN(b_M); |
349 |
cholmod_allocate_dense(cha->ncol, chb->ncol, cha->ncol, chb->xtype, &c); |
CHM_DN chc = cholmod_l_allocate_dense(cha->ncol, chb->ncol, cha->ncol, |
350 |
double alpha[] = {1,0}, beta[] = {0,0}; |
chb->xtype, &c); |
351 |
|
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
352 |
|
double one[] = {1,0}, zero[] = {0,0}; |
353 |
|
R_CheckStack(); |
354 |
|
|
355 |
cholmod_sdmult(cha, 1, alpha, beta, chb, chc, &c); |
cholmod_l_sdmult(cha, 1, one, zero, chb, chc, &c); |
356 |
Free(cha); Free(chb); |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
357 |
UNPROTECT(1); |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1))); |
358 |
return chm_dense_to_SEXP(chc, 1, 0); |
SET_VECTOR_ELT(dn, 1, |
359 |
|
duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1))); |
360 |
|
UNPROTECT(2); |
361 |
|
return chm_dense_to_SEXP(chc, 1, 0, dn); |
362 |
} |
} |
363 |
|
|
364 |
|
/* Computes x'x or x x' -- *also* for Tsparse (triplet = TRUE) |
365 |
|
see Csparse_Csparse_crossprod above for x'y and x y' */ |
366 |
SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet) |
SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet) |
367 |
{ |
{ |
368 |
int trip = asLogical(triplet), |
int trip = asLogical(triplet), |
369 |
tr = asLogical(trans); /* gets reversed because _aat is tcrossprod */ |
tr = asLogical(trans); /* gets reversed because _aat is tcrossprod */ |
370 |
cholmod_triplet |
#ifdef AS_CHM_DIAGU2N_FIXED_FINALLY |
371 |
*cht = trip ? as_cholmod_triplet(x) : (cholmod_triplet*) NULL; |
CHM_TR cht = trip ? AS_CHM_TR(x) : (CHM_TR) NULL; |
372 |
cholmod_sparse *chcp, *chxt, |
#else /* workaround needed:*/ |
373 |
*chx = trip ? cholmod_triplet_to_sparse(cht, cht->nnz, &c) |
CHM_TR cht = trip ? AS_CHM_TR__(Tsparse_diagU2N(x)) : (CHM_TR) NULL; |
374 |
: as_cholmod_sparse(x); |
#endif |
375 |
|
CHM_SP chcp, chxt, |
376 |
|
chx = (trip ? |
377 |
|
cholmod_l_triplet_to_sparse(cht, cht->nnz, &c) : |
378 |
|
AS_CHM_SP(x)); |
379 |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
380 |
|
R_CheckStack(); |
381 |
|
|
382 |
if (!tr) |
if (!tr) chxt = cholmod_l_transpose(chx, chx->xtype, &c); |
383 |
chxt = cholmod_transpose(chx, chx->xtype, &c); |
chcp = cholmod_l_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c); |
384 |
chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c); |
if(!chcp) { |
385 |
if(!chcp) |
UNPROTECT(1); |
386 |
error(_("Csparse_crossprod(): error return from cholmod_aat()")); |
error(_("Csparse_crossprod(): error return from cholmod_l_aat()")); |
|
cholmod_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c); |
|
|
chcp->stype = 1; |
|
|
if (trip) { |
|
|
cholmod_free_sparse(&chx, &c); |
|
|
Free(cht); |
|
|
} else { |
|
|
Free(chx); |
|
387 |
} |
} |
388 |
if (!tr) cholmod_free_sparse(&chxt, &c); |
cholmod_l_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c); |
389 |
/* create dimnames */ |
chcp->stype = 1; |
390 |
SET_VECTOR_ELT(dn, 0, |
if (trip) cholmod_l_free_sparse(&chx, &c); |
391 |
|
if (!tr) cholmod_l_free_sparse(&chxt, &c); |
392 |
|
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
393 |
duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), |
duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), |
394 |
(tr) ? 1 : 0))); |
(tr) ? 0 : 1))); |
395 |
SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(dn, 0))); |
SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(dn, 0))); |
396 |
UNPROTECT(1); |
UNPROTECT(1); |
397 |
return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn); |
return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn); |
399 |
|
|
400 |
SEXP Csparse_drop(SEXP x, SEXP tol) |
SEXP Csparse_drop(SEXP x, SEXP tol) |
401 |
{ |
{ |
402 |
cholmod_sparse *chx = as_cholmod_sparse(x), |
const char *cl = class_P(x); |
403 |
*ans = cholmod_copy(chx, chx->stype, chx->xtype, &c); |
/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */ |
404 |
|
int tr = (cl[1] == 't'); |
405 |
|
CHM_SP chx = AS_CHM_SP__(x); |
406 |
|
CHM_SP ans = cholmod_l_copy(chx, chx->stype, chx->xtype, &c); |
407 |
double dtol = asReal(tol); |
double dtol = asReal(tol); |
408 |
int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
409 |
|
R_CheckStack(); |
410 |
|
|
411 |
if(!cholmod_drop(dtol, ans, &c)) |
if(!cholmod_l_drop(dtol, ans, &c)) |
412 |
error(_("cholmod_drop() failed")); |
error(_("cholmod_l_drop() failed")); |
413 |
Free(chx); |
return chm_sparse_to_SEXP(ans, 1, |
414 |
/* FIXME: currently drops dimnames */ |
tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, |
415 |
return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue); |
Rkind, tr ? diag_P(x) : "", |
416 |
|
GET_SLOT(x, Matrix_DimNamesSym)); |
417 |
} |
} |
418 |
|
|
|
|
|
419 |
SEXP Csparse_horzcat(SEXP x, SEXP y) |
SEXP Csparse_horzcat(SEXP x, SEXP y) |
420 |
{ |
{ |
421 |
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; |
|
422 |
int Rkind = 0; /* only for "d" - FIXME */ |
int Rkind = 0; /* only for "d" - FIXME */ |
423 |
|
R_CheckStack(); |
424 |
|
|
|
ans = cholmod_horzcat(chx, chy, 1, &c); |
|
|
Free(chx); Free(chy); |
|
425 |
/* FIXME: currently drops dimnames */ |
/* FIXME: currently drops dimnames */ |
426 |
return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue); |
return chm_sparse_to_SEXP(cholmod_l_horzcat(chx, chy, 1, &c), |
427 |
|
1, 0, Rkind, "", R_NilValue); |
428 |
} |
} |
429 |
|
|
430 |
SEXP Csparse_vertcat(SEXP x, SEXP y) |
SEXP Csparse_vertcat(SEXP x, SEXP y) |
431 |
{ |
{ |
432 |
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; |
|
433 |
int Rkind = 0; /* only for "d" - FIXME */ |
int Rkind = 0; /* only for "d" - FIXME */ |
434 |
|
R_CheckStack(); |
435 |
|
|
|
ans = cholmod_vertcat(chx, chy, 1, &c); |
|
|
Free(chx); Free(chy); |
|
436 |
/* FIXME: currently drops dimnames */ |
/* FIXME: currently drops dimnames */ |
437 |
return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue); |
return chm_sparse_to_SEXP(cholmod_l_vertcat(chx, chy, 1, &c), |
438 |
|
1, 0, Rkind, "", R_NilValue); |
439 |
} |
} |
440 |
|
|
441 |
SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2) |
SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2) |
442 |
{ |
{ |
443 |
cholmod_sparse *chx = as_cholmod_sparse(x), *ans; |
CHM_SP chx = AS_CHM_SP__(x); |
444 |
int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
445 |
|
CHM_SP ans = cholmod_l_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c); |
446 |
|
R_CheckStack(); |
447 |
|
|
448 |
ans = cholmod_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c); |
return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", |
449 |
Free(chx); |
GET_SLOT(x, Matrix_DimNamesSym)); |
|
return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue); |
|
450 |
} |
} |
451 |
|
|
452 |
SEXP Csparse_diagU2N(SEXP x) |
SEXP Csparse_diagU2N(SEXP x) |
453 |
{ |
{ |
454 |
cholmod_sparse *chx = as_cholmod_sparse(x); |
const char *cl = class_P(x); |
455 |
cholmod_sparse *eye = cholmod_speye(chx->nrow, chx->ncol, chx->xtype, &c); |
/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */ |
456 |
|
if (cl[1] != 't' || *diag_P(x) != 'U') { |
457 |
|
/* "trivially fast" when not triangular (<==> no 'diag' slot), |
458 |
|
or not *unit* triangular */ |
459 |
|
return (x); |
460 |
|
} |
461 |
|
else { /* unit triangular (diag='U'): "fill the diagonal" & diag:= "N" */ |
462 |
|
CHM_SP chx = AS_CHM_SP__(x); |
463 |
|
CHM_SP eye = cholmod_l_speye(chx->nrow, chx->ncol, chx->xtype, &c); |
464 |
double one[] = {1, 0}; |
double one[] = {1, 0}; |
465 |
cholmod_sparse *ans = cholmod_add(chx, eye, one, one, TRUE, TRUE, &c); |
CHM_SP ans = cholmod_l_add(chx, eye, one, one, TRUE, TRUE, &c); |
466 |
int uploT = (strcmp(CHAR(asChar(GET_SLOT(x, Matrix_uploSym))), "U")) ? |
int uploT = (*uplo_P(x) == 'U') ? 1 : -1; |
467 |
-1 : 1; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
|
int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
|
468 |
|
|
469 |
Free(chx); cholmod_free_sparse(&eye, &c); |
R_CheckStack(); |
470 |
|
cholmod_l_free_sparse(&eye, &c); |
471 |
return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N", |
return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N", |
472 |
duplicate(GET_SLOT(x, Matrix_DimNamesSym))); |
GET_SLOT(x, Matrix_DimNamesSym)); |
473 |
|
} |
474 |
|
} |
475 |
|
|
476 |
|
SEXP Csparse_diagN2U(SEXP x) |
477 |
|
{ |
478 |
|
const char *cl = class_P(x); |
479 |
|
/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */ |
480 |
|
if (cl[1] != 't' || *diag_P(x) != 'N') { |
481 |
|
/* "trivially fast" when not triangular (<==> no 'diag' slot), |
482 |
|
or already *unit* triangular */ |
483 |
|
return (x); |
484 |
|
} |
485 |
|
else { /* triangular with diag='N'): now drop the diagonal */ |
486 |
|
/* duplicate, since chx will be modified: */ |
487 |
|
CHM_SP chx = AS_CHM_SP__(duplicate(x)); |
488 |
|
int uploT = (*uplo_P(x) == 'U') ? 1 : -1, |
489 |
|
Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
490 |
|
R_CheckStack(); |
491 |
|
|
492 |
|
chm_diagN2U(chx, uploT, /* do_realloc */ FALSE); |
493 |
|
|
494 |
|
return chm_sparse_to_SEXP(chx, /*dofree*/ 0/* or 1 ?? */, |
495 |
|
uploT, Rkind, "U", |
496 |
|
GET_SLOT(x, Matrix_DimNamesSym)); |
497 |
|
} |
498 |
} |
} |
499 |
|
|
500 |
SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j) |
SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j) |
501 |
{ |
{ |
502 |
cholmod_sparse *chx = as_cholmod_sparse(x); |
CHM_SP chx = AS_CHM_SP__(x); |
503 |
int rsize = (isNull(i)) ? -1 : LENGTH(i), |
int rsize = (isNull(i)) ? -1 : LENGTH(i), |
504 |
csize = (isNull(j)) ? -1 : LENGTH(j); |
csize = (isNull(j)) ? -1 : LENGTH(j); |
505 |
int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
506 |
|
R_CheckStack(); |
507 |
|
|
508 |
if (rsize >= 0 && !isInteger(i)) |
if (rsize >= 0 && !isInteger(i)) |
509 |
error(_("Index i must be NULL or integer")); |
error(_("Index i must be NULL or integer")); |
510 |
if (csize >= 0 && !isInteger(j)) |
if (csize >= 0 && !isInteger(j)) |
511 |
error(_("Index j must be NULL or integer")); |
error(_("Index j must be NULL or integer")); |
512 |
return chm_sparse_to_SEXP(cholmod_submatrix(chx, INTEGER(i), rsize, |
|
513 |
|
return chm_sparse_to_SEXP(cholmod_l_submatrix(chx, INTEGER(i), rsize, |
514 |
INTEGER(j), csize, |
INTEGER(j), csize, |
515 |
TRUE, TRUE, &c), |
TRUE, TRUE, &c), |
516 |
1, 0, Rkind, "", R_NilValue); |
1, 0, Rkind, "", |
517 |
|
/* FIXME: drops dimnames */ R_NilValue); |
518 |
|
} |
519 |
|
|
520 |
|
SEXP Csparse_MatrixMarket(SEXP x, SEXP fname) |
521 |
|
{ |
522 |
|
FILE *f = fopen(CHAR(asChar(fname)), "w"); |
523 |
|
|
524 |
|
if (!f) |
525 |
|
error(_("failure to open file \"%s\" for writing"), |
526 |
|
CHAR(asChar(fname))); |
527 |
|
if (!cholmod_l_write_sparse(f, AS_CHM_SP(x), |
528 |
|
(CHM_SP)NULL, (char*) NULL, &c)) |
529 |
|
error(_("cholmod_l_write_sparse returned error code")); |
530 |
|
fclose(f); |
531 |
|
return R_NilValue; |
532 |
|
} |
533 |
|
|
534 |
|
|
535 |
|
/** |
536 |
|
* Extract the diagonal entries from *triangular* Csparse matrix __or__ a |
537 |
|
* cholmod_sparse factor (LDL = TRUE). |
538 |
|
* |
539 |
|
* @param n dimension of the matrix. |
540 |
|
* @param x_p 'p' (column pointer) slot contents |
541 |
|
* @param x_x 'x' (non-zero entries) slot contents |
542 |
|
* @param perm 'perm' (= permutation vector) slot contents; only used for "diagBack" |
543 |
|
* @param resultKind a (SEXP) string indicating which kind of result is desired. |
544 |
|
* |
545 |
|
* @return a SEXP, either a (double) number or a length n-vector of diagonal entries |
546 |
|
*/ |
547 |
|
SEXP diag_tC_ptr(int n, int *x_p, double *x_x, int *perm, SEXP resultKind) |
548 |
|
/* ^^^^^^ FIXME[Generalize] to int / ... */ |
549 |
|
{ |
550 |
|
const char* res_ch = CHAR(STRING_ELT(resultKind,0)); |
551 |
|
enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log |
552 |
|
} res_kind = ((!strcmp(res_ch, "trace")) ? trace : |
553 |
|
((!strcmp(res_ch, "sumLog")) ? sum_log : |
554 |
|
((!strcmp(res_ch, "prod")) ? prod : |
555 |
|
((!strcmp(res_ch, "diag")) ? diag : |
556 |
|
((!strcmp(res_ch, "diagBack")) ? diag_backpermuted : |
557 |
|
-1))))); |
558 |
|
int i, n_x, i_from = 0; |
559 |
|
SEXP ans = PROTECT(allocVector(REALSXP, |
560 |
|
/* ^^^^ FIXME[Generalize] */ |
561 |
|
(res_kind == diag || |
562 |
|
res_kind == diag_backpermuted) ? n : 1)); |
563 |
|
double *v = REAL(ans); |
564 |
|
/* ^^^^^^ ^^^^ FIXME[Generalize] */ |
565 |
|
|
566 |
|
#define for_DIAG(v_ASSIGN) \ |
567 |
|
for(i = 0; i < n; i++, i_from += n_x) { \ |
568 |
|
/* looking at i-th column */ \ |
569 |
|
n_x = x_p[i+1] - x_p[i];/* #{entries} in this column */ \ |
570 |
|
v_ASSIGN; \ |
571 |
|
} |
572 |
|
|
573 |
|
/* NOTA BENE: we assume -- uplo = "L" i.e. lower triangular matrix |
574 |
|
* for uplo = "U" (makes sense with a "dtCMatrix" !), |
575 |
|
* should use x_x[i_from + (nx - 1)] instead of x_x[i_from], |
576 |
|
* where nx = (x_p[i+1] - x_p[i]) |
577 |
|
*/ |
578 |
|
|
579 |
|
switch(res_kind) { |
580 |
|
case trace: |
581 |
|
v[0] = 0.; |
582 |
|
for_DIAG(v[0] += x_x[i_from]); |
583 |
|
break; |
584 |
|
|
585 |
|
case sum_log: |
586 |
|
v[0] = 0.; |
587 |
|
for_DIAG(v[0] += log(x_x[i_from])); |
588 |
|
break; |
589 |
|
|
590 |
|
case prod: |
591 |
|
v[0] = 1.; |
592 |
|
for_DIAG(v[0] *= x_x[i_from]); |
593 |
|
break; |
594 |
|
|
595 |
|
case diag: |
596 |
|
for_DIAG(v[i] = x_x[i_from]); |
597 |
|
break; |
598 |
|
|
599 |
|
case diag_backpermuted: |
600 |
|
for_DIAG(v[i] = x_x[i_from]); |
601 |
|
|
602 |
|
warning(_("resultKind = 'diagBack' (back-permuted) is experimental")); |
603 |
|
/* now back_permute : */ |
604 |
|
for(i = 0; i < n; i++) { |
605 |
|
double tmp = v[i]; v[i] = v[perm[i]]; v[perm[i]] = tmp; |
606 |
|
/*^^^^ FIXME[Generalize] */ |
607 |
|
} |
608 |
|
break; |
609 |
|
|
610 |
|
default: /* -1 from above */ |
611 |
|
error(_("diag_tC(): invalid 'resultKind'")); |
612 |
|
/* Wall: */ ans = R_NilValue; v = REAL(ans); |
613 |
|
} |
614 |
|
|
615 |
|
UNPROTECT(1); |
616 |
|
return ans; |
617 |
|
} |
618 |
|
|
619 |
|
/** |
620 |
|
* Extract the diagonal entries from *triangular* Csparse matrix __or__ a |
621 |
|
* cholmod_sparse factor (LDL = TRUE). |
622 |
|
* |
623 |
|
* @param pslot 'p' (column pointer) slot of Csparse matrix/factor |
624 |
|
* @param xslot 'x' (non-zero entries) slot of Csparse matrix/factor |
625 |
|
* @param perm_slot 'perm' (= permutation vector) slot of corresponding CHMfactor; |
626 |
|
* only used for "diagBack" |
627 |
|
* @param resultKind a (SEXP) string indicating which kind of result is desired. |
628 |
|
* |
629 |
|
* @return a SEXP, either a (double) number or a length n-vector of diagonal entries |
630 |
|
*/ |
631 |
|
SEXP diag_tC(SEXP pslot, SEXP xslot, SEXP perm_slot, SEXP resultKind) |
632 |
|
{ |
633 |
|
int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */ |
634 |
|
*x_p = INTEGER(pslot), |
635 |
|
*perm = INTEGER(perm_slot); |
636 |
|
double *x_x = REAL(xslot); |
637 |
|
/* ^^^^^^ ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/ |
638 |
|
|
639 |
|
return diag_tC_ptr(n, x_p, x_x, perm, resultKind); |
640 |
|
} |
641 |
|
|
642 |
|
/** |
643 |
|
* Create a Csparse matrix object from indices and/or pointers. |
644 |
|
* |
645 |
|
* @param cls name of actual class of object to create |
646 |
|
* @param i optional integer vector of length nnz of row indices |
647 |
|
* @param j optional integer vector of length nnz of column indices |
648 |
|
* @param p optional integer vector of length np of row or column pointers |
649 |
|
* @param np length of integer vector p. Must be zero if p == (int*)NULL |
650 |
|
* @param x optional vector of values |
651 |
|
* @param nnz length of vectors i, j and/or x, whichever is to be used |
652 |
|
* @param dims optional integer vector of length 2 to be used as |
653 |
|
* dimensions. If dims == (int*)NULL then the maximum row and column |
654 |
|
* index are used as the dimensions. |
655 |
|
* @param dimnames optional list of length 2 to be used as dimnames |
656 |
|
* @param index1 indicator of 1-based indices |
657 |
|
* |
658 |
|
* @return an SEXP of class cls inheriting from CsparseMatrix. |
659 |
|
*/ |
660 |
|
SEXP create_Csparse(char* cls, int* i, int* j, int* p, int np, |
661 |
|
void* x, int nnz, int* dims, SEXP dimnames, |
662 |
|
int index1) |
663 |
|
{ |
664 |
|
SEXP ans; |
665 |
|
int *ij = (int*)NULL, *tri, *trj, |
666 |
|
mi, mj, mp, nrow = -1, ncol = -1; |
667 |
|
int xtype = -1; /* -Wall */ |
668 |
|
CHM_TR T; |
669 |
|
CHM_SP A; |
670 |
|
|
671 |
|
if (np < 0 || nnz < 0) |
672 |
|
error(_("negative vector lengths not allowed: np = %d, nnz = %d"), |
673 |
|
np, nnz); |
674 |
|
if (1 != ((mi = (i == (int*)NULL)) + |
675 |
|
(mj = (j == (int*)NULL)) + |
676 |
|
(mp = (p == (int*)NULL)))) |
677 |
|
error(_("exactly 1 of 'i', 'j' or 'p' must be NULL")); |
678 |
|
if (mp) { |
679 |
|
if (np) error(_("np = %d, must be zero when p is NULL"), np); |
680 |
|
} else { |
681 |
|
if (np) { /* Expand p to form i or j */ |
682 |
|
if (!(p[0])) error(_("p[0] = %d, should be zero"), p[0]); |
683 |
|
for (int ii = 0; ii < np; ii++) |
684 |
|
if (p[ii] > p[ii + 1]) |
685 |
|
error(_("p must be non-decreasing")); |
686 |
|
if (p[np] != nnz) |
687 |
|
error("p[np] = %d != nnz = %d", p[np], nnz); |
688 |
|
ij = Calloc(nnz, int); |
689 |
|
if (mi) { |
690 |
|
i = ij; |
691 |
|
nrow = np; |
692 |
|
} else { |
693 |
|
j = ij; |
694 |
|
ncol = np; |
695 |
|
} |
696 |
|
/* Expand p to 0-based indices */ |
697 |
|
for (int ii = 0; ii < np; ii++) |
698 |
|
for (int jj = p[ii]; jj < p[ii + 1]; jj++) ij[jj] = ii; |
699 |
|
} else { |
700 |
|
if (nnz) |
701 |
|
error(_("Inconsistent dimensions: np = 0 and nnz = %d"), |
702 |
|
nnz); |
703 |
|
} |
704 |
|
} |
705 |
|
/* calculate nrow and ncol */ |
706 |
|
if (nrow < 0) { |
707 |
|
for (int ii = 0; ii < nnz; ii++) { |
708 |
|
int i1 = i[ii] + (index1 ? 0 : 1); /* 1-based index */ |
709 |
|
if (i1 < 1) error(_("invalid row index at position %d"), ii); |
710 |
|
if (i1 > nrow) nrow = i1; |
711 |
|
} |
712 |
|
} |
713 |
|
if (ncol < 0) { |
714 |
|
for (int jj = 0; jj < nnz; jj++) { |
715 |
|
int j1 = j[jj] + (index1 ? 0 : 1); |
716 |
|
if (j1 < 1) error(_("invalid column index at position %d"), jj); |
717 |
|
if (j1 > ncol) ncol = j1; |
718 |
|
} |
719 |
|
} |
720 |
|
if (dims != (int*)NULL) { |
721 |
|
if (dims[0] > nrow) nrow = dims[0]; |
722 |
|
if (dims[1] > ncol) ncol = dims[1]; |
723 |
|
} |
724 |
|
/* check the class name */ |
725 |
|
if (strlen(cls) != 8) |
726 |
|
error(_("strlen of cls argument = %d, should be 8"), strlen(cls)); |
727 |
|
if (!strcmp(cls + 2, "CMatrix")) |
728 |
|
error(_("cls = \"%s\" does not end in \"CMatrix\""), cls); |
729 |
|
switch(cls[0]) { |
730 |
|
case 'd': |
731 |
|
case 'l': |
732 |
|
xtype = CHOLMOD_REAL; |
733 |
|
break; |
734 |
|
case 'n': |
735 |
|
xtype = CHOLMOD_PATTERN; |
736 |
|
break; |
737 |
|
default: |
738 |
|
error(_("cls = \"%s\" must begin with 'd', 'l' or 'n'"), cls); |
739 |
|
} |
740 |
|
if (cls[1] != 'g') |
741 |
|
error(_("Only 'g'eneral sparse matrix types allowed")); |
742 |
|
/* allocate and populate the triplet */ |
743 |
|
T = cholmod_l_allocate_triplet((size_t)nrow, (size_t)ncol, (size_t)nnz, 0, |
744 |
|
xtype, &c); |
745 |
|
T->x = x; |
746 |
|
tri = (int*)T->i; |
747 |
|
trj = (int*)T->j; |
748 |
|
for (int ii = 0; ii < nnz; ii++) { |
749 |
|
tri[ii] = i[ii] - ((!mi && index1) ? 1 : 0); |
750 |
|
trj[ii] = j[ii] - ((!mj && index1) ? 1 : 0); |
751 |
|
} |
752 |
|
/* create the cholmod_sparse structure */ |
753 |
|
A = cholmod_l_triplet_to_sparse(T, nnz, &c); |
754 |
|
cholmod_l_free_triplet(&T, &c); |
755 |
|
/* copy the information to the SEXP */ |
756 |
|
ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cls))); |
757 |
|
/* FIXME: This has been copied from chm_sparse_to_SEXP in chm_common.c */ |
758 |
|
/* allocate and copy common slots */ |
759 |
|
nnz = cholmod_l_nnz(A, &c); |
760 |
|
dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)); |
761 |
|
dims[0] = A->nrow; dims[1] = A->ncol; |
762 |
|
Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, A->ncol + 1)), (int*)A->p, A->ncol + 1); |
763 |
|
Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)), (int*)A->i, nnz); |
764 |
|
switch(cls[1]) { |
765 |
|
case 'd': |
766 |
|
Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)), (double*)A->x, nnz); |
767 |
|
break; |
768 |
|
case 'l': |
769 |
|
error(_("code not yet written for cls = \"lgCMatrix\"")); |
770 |
|
} |
771 |
|
cholmod_l_free_sparse(&A, &c); |
772 |
|
UNPROTECT(1); |
773 |
|
return ans; |
774 |
} |
} |