136 |
} |
} |
137 |
} |
} |
138 |
if (!sorted) |
if (!sorted) |
139 |
/* cannot easily use cholmod_l_sort(.) ... -> "error out" :*/ |
/* cannot easily use cholmod_sort(.) ... -> "error out" :*/ |
140 |
return mkString(_("slot j is not increasing inside a column")); |
return mkString(_("slot j is not increasing inside a column")); |
141 |
else if(!strictly) /* sorted, but not strictly */ |
else if(!strictly) /* sorted, but not strictly */ |
142 |
return mkString(_("slot j is not *strictly* increasing inside a column")); |
return mkString(_("slot j is not *strictly* increasing inside a column")); |
154 |
/* This loses the symmetry property, since cholmod_dense has none, |
/* This loses the symmetry property, since cholmod_dense has none, |
155 |
* BUT, much worse (FIXME!), it also transforms CHOLMOD_PATTERN ("n") matrices |
* BUT, much worse (FIXME!), it also transforms CHOLMOD_PATTERN ("n") matrices |
156 |
* to numeric (CHOLMOD_REAL) ones : */ |
* to numeric (CHOLMOD_REAL) ones : */ |
157 |
CHM_DN chxd = cholmod_l_sparse_to_dense(chxs, &c); |
CHM_DN chxd = cholmod_sparse_to_dense(chxs, &c); |
158 |
int Rkind = (chxs->xtype == CHOLMOD_PATTERN)? -1 : Real_kind(x); |
int Rkind = (chxs->xtype == CHOLMOD_PATTERN)? -1 : Real_kind(x); |
159 |
R_CheckStack(); |
R_CheckStack(); |
160 |
|
|
161 |
return chm_dense_to_SEXP(chxd, 1, Rkind, GET_SLOT(x, Matrix_DimNamesSym)); |
return chm_dense_to_SEXP(chxd, 1, Rkind, GET_SLOT(x, Matrix_DimNamesSym)); |
162 |
} |
} |
163 |
|
|
164 |
|
// FIXME: do not go via CHM (should not be too hard, to just *drop* the x-slot, right? |
165 |
SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri) |
SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri) |
166 |
{ |
{ |
167 |
CHM_SP chxs = AS_CHM_SP__(x); |
CHM_SP chxs = AS_CHM_SP__(x); |
168 |
CHM_SP chxcp = cholmod_l_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c); |
CHM_SP chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c); |
169 |
int tr = asLogical(tri); |
int tr = asLogical(tri); |
170 |
R_CheckStack(); |
R_CheckStack(); |
171 |
|
|
175 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
176 |
} |
} |
177 |
|
|
178 |
|
// n.CMatrix --> [dli].CMatrix (not going through CHM!) |
179 |
|
SEXP nz_pattern_to_Csparse(SEXP x, SEXP res_kind) |
180 |
|
{ |
181 |
|
return nz2Csparse(x, asInteger(res_kind)); |
182 |
|
} |
183 |
|
// n.CMatrix --> [dli].CMatrix (not going through CHM!) |
184 |
|
SEXP nz2Csparse(SEXP x, enum x_slot_kind r_kind) |
185 |
|
{ |
186 |
|
const char *cl_x = class_P(x); |
187 |
|
if(cl_x[0] != 'n') error(_("not a 'n.CMatrix'")); |
188 |
|
if(cl_x[2] != 'C') error(_("not a CsparseMatrix")); |
189 |
|
int nnz = LENGTH(GET_SLOT(x, Matrix_iSym)); |
190 |
|
SEXP ans; |
191 |
|
char *ncl = strdup(cl_x); |
192 |
|
double *dx_x; int *ix_x; |
193 |
|
ncl[0] = (r_kind == x_double ? 'd' : |
194 |
|
(r_kind == x_logical ? 'l' : |
195 |
|
/* else (for now): r_kind == x_integer : */ 'i')); |
196 |
|
PROTECT(ans = NEW_OBJECT(MAKE_CLASS(ncl))); |
197 |
|
// create a correct 'x' slot: |
198 |
|
switch(r_kind) { |
199 |
|
int i; |
200 |
|
case x_double: // 'd' |
201 |
|
dx_x = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)); |
202 |
|
for (i=0; i < nnz; i++) dx_x[i] = 1.; |
203 |
|
break; |
204 |
|
case x_logical: // 'l' |
205 |
|
ix_x = LOGICAL(ALLOC_SLOT(ans, Matrix_xSym, LGLSXP, nnz)); |
206 |
|
for (i=0; i < nnz; i++) ix_x[i] = TRUE; |
207 |
|
break; |
208 |
|
case x_integer: // 'i' |
209 |
|
ix_x = INTEGER(ALLOC_SLOT(ans, Matrix_xSym, INTSXP, nnz)); |
210 |
|
for (i=0; i < nnz; i++) ix_x[i] = 1; |
211 |
|
break; |
212 |
|
|
213 |
|
default: |
214 |
|
error(_("nz2Csparse(): invalid/non-implemented r_kind = %d"), |
215 |
|
r_kind); |
216 |
|
} |
217 |
|
|
218 |
|
// now copy all other slots : |
219 |
|
slot_dup(ans, x, Matrix_iSym); |
220 |
|
slot_dup(ans, x, Matrix_pSym); |
221 |
|
slot_dup(ans, x, Matrix_DimSym); |
222 |
|
slot_dup(ans, x, Matrix_DimNamesSym); |
223 |
|
if(ncl[1] != 'g') { // symmetric or triangular ... |
224 |
|
slot_dup_if_has(ans, x, Matrix_uploSym); |
225 |
|
slot_dup_if_has(ans, x, Matrix_diagSym); |
226 |
|
} |
227 |
|
UNPROTECT(1); |
228 |
|
return ans; |
229 |
|
} |
230 |
|
|
231 |
SEXP Csparse_to_matrix(SEXP x) |
SEXP Csparse_to_matrix(SEXP x) |
232 |
{ |
{ |
233 |
return chm_dense_to_matrix(cholmod_l_sparse_to_dense(AS_CHM_SP__(x), &c), |
return chm_dense_to_matrix(cholmod_sparse_to_dense(AS_CHM_SP__(x), &c), |
234 |
1 /*do_free*/, GET_SLOT(x, Matrix_DimNamesSym)); |
1 /*do_free*/, GET_SLOT(x, Matrix_DimNamesSym)); |
235 |
} |
} |
236 |
|
|
237 |
SEXP Csparse_to_Tsparse(SEXP x, SEXP tri) |
SEXP Csparse_to_Tsparse(SEXP x, SEXP tri) |
238 |
{ |
{ |
239 |
CHM_SP chxs = AS_CHM_SP__(x); |
CHM_SP chxs = AS_CHM_SP__(x); |
240 |
CHM_TR chxt = cholmod_l_sparse_to_triplet(chxs, &c); |
CHM_TR chxt = cholmod_sparse_to_triplet(chxs, &c); |
241 |
int tr = asLogical(tri); |
int tr = asLogical(tri); |
242 |
int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
243 |
R_CheckStack(); |
R_CheckStack(); |
257 |
|
|
258 |
if (!(chx->stype)) |
if (!(chx->stype)) |
259 |
error(_("Nonsymmetric matrix in Csparse_symmetric_to_general")); |
error(_("Nonsymmetric matrix in Csparse_symmetric_to_general")); |
260 |
chgx = cholmod_l_copy(chx, /* stype: */ 0, chx->xtype, &c); |
chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c); |
261 |
/* xtype: pattern, "real", complex or .. */ |
/* xtype: pattern, "real", complex or .. */ |
262 |
return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
263 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
270 |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
271 |
R_CheckStack(); |
R_CheckStack(); |
272 |
|
|
273 |
chgx = cholmod_l_copy(chx, /* stype: */ uploT, chx->xtype, &c); |
chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c); |
274 |
/* xtype: pattern, "real", complex or .. */ |
/* xtype: pattern, "real", complex or .. */ |
275 |
return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
276 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
282 |
* since cholmod (& cs) lacks sparse 'int' matrices */ |
* since cholmod (& cs) lacks sparse 'int' matrices */ |
283 |
CHM_SP chx = AS_CHM_SP__(x); |
CHM_SP chx = AS_CHM_SP__(x); |
284 |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
285 |
CHM_SP chxt = cholmod_l_transpose(chx, chx->xtype, &c); |
CHM_SP chxt = cholmod_transpose(chx, chx->xtype, &c); |
286 |
SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp; |
SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp; |
287 |
int tr = asLogical(tri); |
int tr = asLogical(tri); |
288 |
R_CheckStack(); |
R_CheckStack(); |
301 |
CHM_SP |
CHM_SP |
302 |
cha = AS_CHM_SP(a), |
cha = AS_CHM_SP(a), |
303 |
chb = AS_CHM_SP(b), |
chb = AS_CHM_SP(b), |
304 |
chc = cholmod_l_ssmult(cha, chb, /*out_stype:*/ 0, |
chc = cholmod_ssmult(cha, chb, /*out_stype:*/ 0, |
305 |
/* values:= is_numeric (T/F) */ cha->xtype > 0, |
/* values:= is_numeric (T/F) */ cha->xtype > 0, |
306 |
/*out sorted:*/ 1, &c); |
/*out sorted:*/ 1, &c); |
307 |
const char *cl_a = class_P(a), *cl_b = class_P(b); |
const char *cl_a = class_P(a), *cl_b = class_P(b); |
352 |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
353 |
R_CheckStack(); |
R_CheckStack(); |
354 |
|
|
355 |
chTr = cholmod_l_transpose((tr) ? chb : cha, chb->xtype, &c); |
chTr = cholmod_transpose((tr) ? chb : cha, chb->xtype, &c); |
356 |
chc = cholmod_l_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb, |
chc = cholmod_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb, |
357 |
/*out_stype:*/ 0, cha->xtype, /*out sorted:*/ 1, &c); |
/*out_stype:*/ 0, cha->xtype, /*out sorted:*/ 1, &c); |
358 |
cholmod_l_free_sparse(&chTr, &c); |
cholmod_free_sparse(&chTr, &c); |
359 |
|
|
360 |
/* Preserve triangularity and unit-triangularity if appropriate; |
/* Preserve triangularity and unit-triangularity if appropriate; |
361 |
* see Csparse_Csparse_prod() for comments */ |
* see Csparse_Csparse_prod() for comments */ |
381 |
CHM_SP cha = AS_CHM_SP(a); |
CHM_SP cha = AS_CHM_SP(a); |
382 |
SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
383 |
CHM_DN chb = AS_CHM_DN(b_M); |
CHM_DN chb = AS_CHM_DN(b_M); |
384 |
CHM_DN chc = cholmod_l_allocate_dense(cha->nrow, chb->ncol, cha->nrow, |
CHM_DN chc = cholmod_allocate_dense(cha->nrow, chb->ncol, cha->nrow, |
385 |
chb->xtype, &c); |
chb->xtype, &c); |
386 |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
387 |
double one[] = {1,0}, zero[] = {0,0}; |
double one[] = {1,0}, zero[] = {0,0}; |
388 |
|
int nprot = 2; |
389 |
R_CheckStack(); |
R_CheckStack(); |
390 |
|
/* Tim Davis, please FIXME: currently (2010-11) *fails* when a is a pattern matrix:*/ |
391 |
cholmod_l_sdmult(cha, 0, one, zero, chb, chc, &c); |
if(cha->xtype == CHOLMOD_PATTERN) { |
392 |
|
/* warning(_("Csparse_dense_prod(): cholmod_sdmult() not yet implemented for pattern./ ngCMatrix" */ |
393 |
|
/* " --> slightly inefficient coercion")); */ |
394 |
|
|
395 |
|
// This *fails* to produce a CHOLMOD_REAL .. |
396 |
|
// CHM_SP chd = cholmod_l_copy(cha, cha->stype, CHOLMOD_REAL, &c); |
397 |
|
// --> use our Matrix-classes |
398 |
|
SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++; |
399 |
|
cha = AS_CHM_SP(da); |
400 |
|
} |
401 |
|
cholmod_sdmult(cha, 0, one, zero, chb, chc, &c); |
402 |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
403 |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); |
404 |
SET_VECTOR_ELT(dn, 1, |
SET_VECTOR_ELT(dn, 1, |
405 |
duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1))); |
duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1))); |
406 |
UNPROTECT(2); |
UNPROTECT(nprot); |
407 |
return chm_dense_to_SEXP(chc, 1, 0, dn); |
return chm_dense_to_SEXP(chc, 1, 0, dn); |
408 |
} |
} |
409 |
|
|
412 |
CHM_SP cha = AS_CHM_SP(a); |
CHM_SP cha = AS_CHM_SP(a); |
413 |
SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
414 |
CHM_DN chb = AS_CHM_DN(b_M); |
CHM_DN chb = AS_CHM_DN(b_M); |
415 |
CHM_DN chc = cholmod_l_allocate_dense(cha->ncol, chb->ncol, cha->ncol, |
CHM_DN chc = cholmod_allocate_dense(cha->ncol, chb->ncol, cha->ncol, |
416 |
chb->xtype, &c); |
chb->xtype, &c); |
417 |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); int nprot = 2; |
418 |
double one[] = {1,0}, zero[] = {0,0}; |
double one[] = {1,0}, zero[] = {0,0}; |
419 |
R_CheckStack(); |
R_CheckStack(); |
420 |
|
// -- see Csparse_dense_prod() above : |
421 |
cholmod_l_sdmult(cha, 1, one, zero, chb, chc, &c); |
if(cha->xtype == CHOLMOD_PATTERN) { |
422 |
|
SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++; |
423 |
|
cha = AS_CHM_SP(da); |
424 |
|
} |
425 |
|
cholmod_sdmult(cha, 1, one, zero, chb, chc, &c); |
426 |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
427 |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1))); |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1))); |
428 |
SET_VECTOR_ELT(dn, 1, |
SET_VECTOR_ELT(dn, 1, |
429 |
duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1))); |
duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1))); |
430 |
UNPROTECT(2); |
UNPROTECT(nprot); |
431 |
return chm_dense_to_SEXP(chc, 1, 0, dn); |
return chm_dense_to_SEXP(chc, 1, 0, dn); |
432 |
} |
} |
433 |
|
|
445 |
#endif |
#endif |
446 |
CHM_SP chcp, chxt, |
CHM_SP chcp, chxt, |
447 |
chx = (trip ? |
chx = (trip ? |
448 |
cholmod_l_triplet_to_sparse(cht, cht->nnz, &c) : |
cholmod_triplet_to_sparse(cht, cht->nnz, &c) : |
449 |
AS_CHM_SP(x)); |
AS_CHM_SP(x)); |
450 |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
451 |
R_CheckStack(); |
R_CheckStack(); |
452 |
|
|
453 |
if (!tr) chxt = cholmod_l_transpose(chx, chx->xtype, &c); |
if (!tr) chxt = cholmod_transpose(chx, chx->xtype, &c); |
454 |
chcp = cholmod_l_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c); |
chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c); |
455 |
if(!chcp) { |
if(!chcp) { |
456 |
UNPROTECT(1); |
UNPROTECT(1); |
457 |
error(_("Csparse_crossprod(): error return from cholmod_l_aat()")); |
error(_("Csparse_crossprod(): error return from cholmod_aat()")); |
458 |
} |
} |
459 |
cholmod_l_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c); |
cholmod_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c); |
460 |
chcp->stype = 1; |
chcp->stype = 1; |
461 |
if (trip) cholmod_l_free_sparse(&chx, &c); |
if (trip) cholmod_free_sparse(&chx, &c); |
462 |
if (!tr) cholmod_l_free_sparse(&chxt, &c); |
if (!tr) cholmod_free_sparse(&chxt, &c); |
463 |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
464 |
duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), |
duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), |
465 |
(tr) ? 0 : 1))); |
(tr) ? 0 : 1))); |
472 |
return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn); |
return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn); |
473 |
} |
} |
474 |
|
|
475 |
|
/* Csparse_drop(x, tol): drop entries with absolute value < tol, i.e, |
476 |
|
* at least all "explicit" zeros */ |
477 |
SEXP Csparse_drop(SEXP x, SEXP tol) |
SEXP Csparse_drop(SEXP x, SEXP tol) |
478 |
{ |
{ |
479 |
const char *cl = class_P(x); |
const char *cl = class_P(x); |
480 |
/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */ |
/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */ |
481 |
int tr = (cl[1] == 't'); |
int tr = (cl[1] == 't'); |
482 |
CHM_SP chx = AS_CHM_SP__(x); |
CHM_SP chx = AS_CHM_SP__(x); |
483 |
CHM_SP ans = cholmod_l_copy(chx, chx->stype, chx->xtype, &c); |
CHM_SP ans = cholmod_copy(chx, chx->stype, chx->xtype, &c); |
484 |
double dtol = asReal(tol); |
double dtol = asReal(tol); |
485 |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
486 |
R_CheckStack(); |
R_CheckStack(); |
487 |
|
|
488 |
if(!cholmod_l_drop(dtol, ans, &c)) |
if(!cholmod_drop(dtol, ans, &c)) |
489 |
error(_("cholmod_l_drop() failed")); |
error(_("cholmod_drop() failed")); |
490 |
return chm_sparse_to_SEXP(ans, 1, |
return chm_sparse_to_SEXP(ans, 1, |
491 |
tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, |
tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, |
492 |
Rkind, tr ? diag_P(x) : "", |
Rkind, tr ? diag_P(x) : "", |
502 |
R_CheckStack(); |
R_CheckStack(); |
503 |
|
|
504 |
/* TODO: currently drops dimnames - and we fix at R level */ |
/* TODO: currently drops dimnames - and we fix at R level */ |
505 |
return chm_sparse_to_SEXP(cholmod_l_horzcat(chx, chy, 1, &c), |
return chm_sparse_to_SEXP(cholmod_horzcat(chx, chy, 1, &c), |
506 |
1, 0, Rkind, "", R_NilValue); |
1, 0, Rkind, "", R_NilValue); |
507 |
} |
} |
508 |
|
|
515 |
R_CheckStack(); |
R_CheckStack(); |
516 |
|
|
517 |
/* TODO: currently drops dimnames - and we fix at R level */ |
/* TODO: currently drops dimnames - and we fix at R level */ |
518 |
return chm_sparse_to_SEXP(cholmod_l_vertcat(chx, chy, 1, &c), |
return chm_sparse_to_SEXP(cholmod_vertcat(chx, chy, 1, &c), |
519 |
1, 0, Rkind, "", R_NilValue); |
1, 0, Rkind, "", R_NilValue); |
520 |
} |
} |
521 |
|
|
523 |
{ |
{ |
524 |
CHM_SP chx = AS_CHM_SP__(x); |
CHM_SP chx = AS_CHM_SP__(x); |
525 |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
526 |
CHM_SP ans = cholmod_l_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c); |
CHM_SP ans = cholmod_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c); |
527 |
R_CheckStack(); |
R_CheckStack(); |
528 |
|
|
529 |
return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", |
return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", |
541 |
} |
} |
542 |
else { /* unit triangular (diag='U'): "fill the diagonal" & diag:= "N" */ |
else { /* unit triangular (diag='U'): "fill the diagonal" & diag:= "N" */ |
543 |
CHM_SP chx = AS_CHM_SP__(x); |
CHM_SP chx = AS_CHM_SP__(x); |
544 |
CHM_SP eye = cholmod_l_speye(chx->nrow, chx->ncol, chx->xtype, &c); |
CHM_SP eye = cholmod_speye(chx->nrow, chx->ncol, chx->xtype, &c); |
545 |
double one[] = {1, 0}; |
double one[] = {1, 0}; |
546 |
CHM_SP ans = cholmod_l_add(chx, eye, one, one, TRUE, TRUE, &c); |
CHM_SP ans = cholmod_add(chx, eye, one, one, TRUE, TRUE, &c); |
547 |
int uploT = (*uplo_P(x) == 'U') ? 1 : -1; |
int uploT = (*uplo_P(x) == 'U') ? 1 : -1; |
548 |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
549 |
|
|
550 |
R_CheckStack(); |
R_CheckStack(); |
551 |
cholmod_l_free_sparse(&eye, &c); |
cholmod_free_sparse(&eye, &c); |
552 |
return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N", |
return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N", |
553 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
554 |
} |
} |
602 |
|
|
603 |
if (chx->stype) /* symmetricMatrix */ |
if (chx->stype) /* symmetricMatrix */ |
604 |
/* for now, cholmod_submatrix() only accepts "generalMatrix" */ |
/* for now, cholmod_submatrix() only accepts "generalMatrix" */ |
605 |
chx = cholmod_l_copy(chx, /* stype: */ 0, chx->xtype, &c); |
chx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c); |
606 |
|
|
607 |
return chm_sparse_to_SEXP(cholmod_l_submatrix(chx, |
return chm_sparse_to_SEXP(cholmod_submatrix(chx, |
608 |
(rsize < 0) ? NULL : INTEGER(i), rsize, |
(rsize < 0) ? NULL : INTEGER(i), rsize, |
609 |
(csize < 0) ? NULL : INTEGER(j), csize, |
(csize < 0) ? NULL : INTEGER(j), csize, |
610 |
TRUE, TRUE, &c), |
TRUE, TRUE, &c), |
612 |
/* FIXME: drops dimnames */ R_NilValue); |
/* FIXME: drops dimnames */ R_NilValue); |
613 |
} |
} |
614 |
|
|
615 |
|
/** |
616 |
|
* Subassignment: x[i,j] <- value |
617 |
|
* |
618 |
|
* @param x |
619 |
|
* @param i_ integer row index 0-origin vector (as returned from R .ind.prep2()) |
620 |
|
* @param j_ integer column index 0-origin vector |
621 |
|
* @param value currently must be a dsparseVector {which is recycled if needed} |
622 |
|
* |
623 |
|
* @return a Csparse matrix like x, but with the values replaced |
624 |
|
*/ |
625 |
|
SEXP Csparse_subassign(SEXP x, SEXP i_, SEXP j_, SEXP value) |
626 |
|
{ |
627 |
|
static const char |
628 |
|
*valid_cM [] = {"dgCMatrix",// the only one, for "the moment", more later |
629 |
|
""}, |
630 |
|
*valid_spv[] = {"dsparseVector", |
631 |
|
""}; |
632 |
|
|
633 |
|
int ctype = Matrix_check_class_etc(x, valid_cM); |
634 |
|
if (ctype < 0) |
635 |
|
error(_("invalid class of 'x' in Csparse_subassign()")); |
636 |
|
// value: assume a "dsparseVector" for now -- slots: (i, length, x) |
637 |
|
ctype = Matrix_check_class_etc(value, valid_spv); |
638 |
|
if (ctype < 0) |
639 |
|
error(_("invalid class of 'value' in Csparse_subassign()")); |
640 |
|
|
641 |
|
SEXP ans, |
642 |
|
pslot = GET_SLOT(x, Matrix_pSym), |
643 |
|
islot = GET_SLOT(x, Matrix_iSym), |
644 |
|
i_cp = PROTECT((TYPEOF(i_) == INTSXP) ? |
645 |
|
duplicate(i_) : coerceVector(i_, INTSXP)), |
646 |
|
j_cp = PROTECT((TYPEOF(j_) == INTSXP) ? |
647 |
|
duplicate(j_) : coerceVector(j_, INTSXP)), |
648 |
|
// for d.CMatrix and l.CMatrix but not n.CMatrix: |
649 |
|
xslot = GET_SLOT(x, Matrix_xSym); |
650 |
|
|
651 |
|
int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), |
652 |
|
nrow = dims[0], |
653 |
|
ncol = dims[1], |
654 |
|
*xp = INTEGER(pslot), |
655 |
|
*xi = INTEGER(islot), |
656 |
|
*ii = INTEGER(i_cp), len_i = LENGTH(i_cp), |
657 |
|
*jj = INTEGER(j_cp), len_j = LENGTH(j_cp), |
658 |
|
i, j, k; |
659 |
|
int *val_i = INTEGER(GET_SLOT(value, Matrix_iSym)); |
660 |
|
// for dsparseVector only: |
661 |
|
double *val_x = REAL (GET_SLOT(value, Matrix_xSym)); |
662 |
|
int len_val = asInteger(GET_SLOT(value, Matrix_lengthSym)); |
663 |
|
int p_last = xp[0]; |
664 |
|
|
665 |
|
// for d.CMatrix only: |
666 |
|
double *xx = REAL(xslot); |
667 |
|
double ind; // the index that goes all the way from 1:(len_i * len_j) |
668 |
|
|
669 |
|
PROTECT(ans = duplicate(x)); |
670 |
|
for(j = 0; j < ncol; j++) { |
671 |
|
// FIXME |
672 |
|
// .... |
673 |
|
// .... |
674 |
|
// .... |
675 |
|
// .... |
676 |
|
|
677 |
|
|
678 |
|
|
679 |
|
|
680 |
|
|
681 |
|
|
682 |
|
|
683 |
|
// .... |
684 |
|
// .... |
685 |
|
// .... |
686 |
|
// .... |
687 |
|
// .... |
688 |
|
} |
689 |
|
UNPROTECT(3); |
690 |
|
return ans; |
691 |
|
} |
692 |
|
|
693 |
SEXP Csparse_MatrixMarket(SEXP x, SEXP fname) |
SEXP Csparse_MatrixMarket(SEXP x, SEXP fname) |
694 |
{ |
{ |
695 |
FILE *f = fopen(CHAR(asChar(fname)), "w"); |
FILE *f = fopen(CHAR(asChar(fname)), "w"); |
697 |
if (!f) |
if (!f) |
698 |
error(_("failure to open file \"%s\" for writing"), |
error(_("failure to open file \"%s\" for writing"), |
699 |
CHAR(asChar(fname))); |
CHAR(asChar(fname))); |
700 |
if (!cholmod_l_write_sparse(f, AS_CHM_SP(x), |
if (!cholmod_write_sparse(f, AS_CHM_SP(x), |
701 |
(CHM_SP)NULL, (char*) NULL, &c)) |
(CHM_SP)NULL, (char*) NULL, &c)) |
702 |
error(_("cholmod_l_write_sparse returned error code")); |
error(_("cholmod_write_sparse returned error code")); |
703 |
fclose(f); |
fclose(f); |
704 |
return R_NilValue; |
return R_NilValue; |
705 |
} |
} |
913 |
if (cls[1] != 'g') |
if (cls[1] != 'g') |
914 |
error(_("Only 'g'eneral sparse matrix types allowed")); |
error(_("Only 'g'eneral sparse matrix types allowed")); |
915 |
/* allocate and populate the triplet */ |
/* allocate and populate the triplet */ |
916 |
T = cholmod_l_allocate_triplet((size_t)nrow, (size_t)ncol, (size_t)nnz, 0, |
T = cholmod_allocate_triplet((size_t)nrow, (size_t)ncol, (size_t)nnz, 0, |
917 |
xtype, &c); |
xtype, &c); |
918 |
T->x = x; |
T->x = x; |
919 |
tri = (int*)T->i; |
tri = (int*)T->i; |
923 |
trj[ii] = j[ii] - ((!mj && index1) ? 1 : 0); |
trj[ii] = j[ii] - ((!mj && index1) ? 1 : 0); |
924 |
} |
} |
925 |
/* create the cholmod_sparse structure */ |
/* create the cholmod_sparse structure */ |
926 |
A = cholmod_l_triplet_to_sparse(T, nnz, &c); |
A = cholmod_triplet_to_sparse(T, nnz, &c); |
927 |
cholmod_l_free_triplet(&T, &c); |
cholmod_free_triplet(&T, &c); |
928 |
/* copy the information to the SEXP */ |
/* copy the information to the SEXP */ |
929 |
ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cls))); |
ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cls))); |
930 |
/* FIXME: This has been copied from chm_sparse_to_SEXP in chm_common.c */ |
/* FIXME: This has been copied from chm_sparse_to_SEXP in chm_common.c */ |
931 |
/* allocate and copy common slots */ |
/* allocate and copy common slots */ |
932 |
nnz = cholmod_l_nnz(A, &c); |
nnz = cholmod_nnz(A, &c); |
933 |
dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)); |
dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)); |
934 |
dims[0] = A->nrow; dims[1] = A->ncol; |
dims[0] = A->nrow; dims[1] = A->ncol; |
935 |
Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, A->ncol + 1)), (int*)A->p, A->ncol + 1); |
Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, A->ncol + 1)), (int*)A->p, A->ncol + 1); |
941 |
case 'l': |
case 'l': |
942 |
error(_("code not yet written for cls = \"lgCMatrix\"")); |
error(_("code not yet written for cls = \"lgCMatrix\"")); |
943 |
} |
} |
944 |
cholmod_l_free_sparse(&A, &c); |
/* FIXME: dimnames are *NOT* put there yet (if non-NULL) */ |
945 |
|
cholmod_free_sparse(&A, &c); |
946 |
UNPROTECT(1); |
UNPROTECT(1); |
947 |
return ans; |
return ans; |
948 |
} |
} |