206 |
CHM_SP |
CHM_SP |
207 |
cha = AS_CHM_SP(Csparse_diagU2N(a)), |
cha = AS_CHM_SP(Csparse_diagU2N(a)), |
208 |
chb = AS_CHM_SP(Csparse_diagU2N(b)), |
chb = AS_CHM_SP(Csparse_diagU2N(b)), |
209 |
chc = cholmod_ssmult(cha, chb, 0, cha->xtype, 1, &c); |
chc = cholmod_ssmult(cha, chb, /*out_stype:*/ 0, |
210 |
|
cha->xtype, /*out sorted:*/ 1, &c); |
211 |
|
const char *cl_a = class_P(a), *cl_b = class_P(b); |
212 |
|
char diag[] = {'\0', '\0'}; |
213 |
|
int uploT = 0; |
214 |
SEXP dn = allocVector(VECSXP, 2); |
SEXP dn = allocVector(VECSXP, 2); |
215 |
R_CheckStack(); |
R_CheckStack(); |
216 |
|
|
217 |
|
/* Preserve triangularity and even unit-triangularity if appropriate. |
218 |
|
* Note that in that case, the multiplication itself should happen |
219 |
|
* faster. But there's no support for that in CHOLMOD */ |
220 |
|
|
221 |
|
/* UGLY hack -- rather should have (fast!) C-level version of |
222 |
|
* is(a, "triangularMatrix") etc */ |
223 |
|
if (cl_a[1] == 't' && cl_b[1] == 't') |
224 |
|
/* FIXME: fails for "Cholesky","BunchKaufmann"..*/ |
225 |
|
if(*uplo_P(a) == *uplo_P(b)) { /* both upper, or both lower tri. */ |
226 |
|
uploT = (*uplo_P(a) == 'U') ? 1 : -1; |
227 |
|
if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */ |
228 |
|
/* "remove the diagonal entries": */ |
229 |
|
chm_diagN2U(chc, uploT, /* do_realloc */ FALSE); |
230 |
|
diag[0]= 'U'; |
231 |
|
} |
232 |
|
else diag[0]= 'N'; |
233 |
|
} |
234 |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
235 |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); |
236 |
SET_VECTOR_ELT(dn, 1, |
SET_VECTOR_ELT(dn, 1, |
237 |
duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1))); |
duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1))); |
238 |
return chm_sparse_to_SEXP(chc, 1, 0, 0, "", dn); |
return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn); |
239 |
} |
} |
240 |
|
|
241 |
SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b, SEXP trans) |
SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b, SEXP trans) |
245 |
cha = AS_CHM_SP(Csparse_diagU2N(a)), |
cha = AS_CHM_SP(Csparse_diagU2N(a)), |
246 |
chb = AS_CHM_SP(Csparse_diagU2N(b)), |
chb = AS_CHM_SP(Csparse_diagU2N(b)), |
247 |
chTr, chc; |
chTr, chc; |
248 |
|
const char *cl_a = class_P(a), *cl_b = class_P(b); |
249 |
|
char diag[] = {'\0', '\0'}; |
250 |
|
int uploT = 0; |
251 |
SEXP dn = allocVector(VECSXP, 2); |
SEXP dn = allocVector(VECSXP, 2); |
252 |
R_CheckStack(); |
R_CheckStack(); |
253 |
|
|
254 |
chTr = cholmod_transpose((tr) ? chb : cha, chb->xtype, &c); |
chTr = cholmod_transpose((tr) ? chb : cha, chb->xtype, &c); |
255 |
chc = cholmod_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb, |
chc = cholmod_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb, |
256 |
0, cha->xtype, 1, &c); |
/*out_stype:*/ 0, cha->xtype, /*out sorted:*/ 1, &c); |
257 |
cholmod_free_sparse(&chTr, &c); |
cholmod_free_sparse(&chTr, &c); |
258 |
|
|
259 |
|
/* Preserve triangularity and unit-triangularity if appropriate; |
260 |
|
* see Csparse_Csparse_prod() for comments */ |
261 |
|
if (cl_a[1] == 't' && cl_b[1] == 't') |
262 |
|
if(*uplo_P(a) != *uplo_P(b)) { /* one 'U', the other 'L' */ |
263 |
|
uploT = (*uplo_P(b) == 'U') ? 1 : -1; |
264 |
|
if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */ |
265 |
|
chm_diagN2U(chc, uploT, /* do_realloc */ FALSE); |
266 |
|
diag[0]= 'U'; |
267 |
|
} |
268 |
|
else diag[0]= 'N'; |
269 |
|
} |
270 |
|
|
271 |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
272 |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), (tr) ? 0 : 1))); |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), (tr) ? 0 : 1))); |
273 |
SET_VECTOR_ELT(dn, 1, |
SET_VECTOR_ELT(dn, 1, |
274 |
duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), (tr) ? 0 : 1))); |
duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), (tr) ? 0 : 1))); |
275 |
return chm_sparse_to_SEXP(chc, 1, 0, 0, "", dn); |
return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn); |
276 |
} |
} |
277 |
|
|
278 |
SEXP Csparse_dense_prod(SEXP a, SEXP b) |
SEXP Csparse_dense_prod(SEXP a, SEXP b) |
315 |
return chm_dense_to_SEXP(chc, 1, 0, dn); |
return chm_dense_to_SEXP(chc, 1, 0, dn); |
316 |
} |
} |
317 |
|
|
318 |
/* Computes x'x or x x' -- see Csparse_Csparse_crossprod above for x'y and x y' */ |
/* Computes x'x or x x' -- *also* for Tsparse (triplet = TRUE) |
319 |
|
see Csparse_Csparse_crossprod above for x'y and x y' */ |
320 |
SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet) |
SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet) |
321 |
{ |
{ |
322 |
int trip = asLogical(triplet), |
int trip = asLogical(triplet), |
399 |
const char *cl = class_P(x); |
const char *cl = class_P(x); |
400 |
/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */ |
/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */ |
401 |
if (cl[1] != 't' || *diag_P(x) != 'U') { |
if (cl[1] != 't' || *diag_P(x) != 'U') { |
402 |
/* "trivially fast" when not triangular (<==> no 'diag' slot), or not *unit* triangular */ |
/* "trivially fast" when not triangular (<==> no 'diag' slot), |
403 |
|
or not *unit* triangular */ |
404 |
return (x); |
return (x); |
405 |
} |
} |
406 |
else { |
else { /* unit triangular (diag='U'): "fill the diagonal" & diag:= "N" */ |
407 |
CHM_SP chx = AS_CHM_SP(x); |
CHM_SP chx = AS_CHM_SP(x); |
408 |
CHM_SP eye = cholmod_speye(chx->nrow, chx->ncol, chx->xtype, &c); |
CHM_SP eye = cholmod_speye(chx->nrow, chx->ncol, chx->xtype, &c); |
409 |
double one[] = {1, 0}; |
double one[] = {1, 0}; |
417 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
418 |
} |
} |
419 |
} |
} |
420 |
|
|
421 |
|
SEXP Csparse_diagN2U(SEXP x) |
422 |
|
{ |
423 |
|
const char *cl = class_P(x); |
424 |
|
/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */ |
425 |
|
if (cl[1] != 't' || *diag_P(x) != 'N') { |
426 |
|
/* "trivially fast" when not triangular (<==> no 'diag' slot), |
427 |
|
or already *unit* triangular */ |
428 |
|
return (x); |
429 |
|
} |
430 |
|
else { /* triangular with diag='N'): now drop the diagonal */ |
431 |
|
/* duplicate, since chx will be modified: */ |
432 |
|
CHM_SP chx = AS_CHM_SP(duplicate(x)); |
433 |
|
int uploT = (*uplo_P(x) == 'U') ? 1 : -1, |
434 |
|
Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
435 |
|
R_CheckStack(); |
436 |
|
|
437 |
|
chm_diagN2U(chx, uploT, /* do_realloc */ FALSE); |
438 |
|
|
439 |
|
return chm_sparse_to_SEXP(chx, /*dofree*/ 0/* or 1 ?? */, |
440 |
|
uploT, Rkind, "U", |
441 |
|
GET_SLOT(x, Matrix_DimNamesSym)); |
442 |
|
} |
443 |
|
} |
444 |
|
|
445 |
SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j) |
SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j) |
446 |
{ |
{ |