1 |
/* Sparse matrices in compressed column-oriented form */ |
/* Sparse matrices in compressed column-oriented form */ |
2 |
|
|
3 |
#include "Csparse.h" |
#include "Csparse.h" |
4 |
#include "Tsparse.h" |
#include "Tsparse.h" |
5 |
#include "chm_common.h" |
#include "chm_common.h" |
137 |
} |
} |
138 |
} |
} |
139 |
if (!sorted) |
if (!sorted) |
140 |
/* cannot easily use cholmod_l_sort(.) ... -> "error out" :*/ |
/* cannot easily use cholmod_sort(.) ... -> "error out" :*/ |
141 |
return mkString(_("slot j is not increasing inside a column")); |
return mkString(_("slot j is not increasing inside a column")); |
142 |
else if(!strictly) /* sorted, but not strictly */ |
else if(!strictly) /* sorted, but not strictly */ |
143 |
return mkString(_("slot j is not *strictly* increasing inside a column")); |
return mkString(_("slot j is not *strictly* increasing inside a column")); |
155 |
/* This loses the symmetry property, since cholmod_dense has none, |
/* This loses the symmetry property, since cholmod_dense has none, |
156 |
* BUT, much worse (FIXME!), it also transforms CHOLMOD_PATTERN ("n") matrices |
* BUT, much worse (FIXME!), it also transforms CHOLMOD_PATTERN ("n") matrices |
157 |
* to numeric (CHOLMOD_REAL) ones : */ |
* to numeric (CHOLMOD_REAL) ones : */ |
158 |
CHM_DN chxd = cholmod_l_sparse_to_dense(chxs, &c); |
CHM_DN chxd = cholmod_sparse_to_dense(chxs, &c); |
159 |
int Rkind = (chxs->xtype == CHOLMOD_PATTERN)? -1 : Real_kind(x); |
int Rkind = (chxs->xtype == CHOLMOD_PATTERN)? -1 : Real_kind(x); |
160 |
R_CheckStack(); |
R_CheckStack(); |
161 |
|
|
162 |
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)); |
163 |
} |
} |
164 |
|
|
165 |
|
// FIXME: do not go via CHM (should not be too hard, to just *drop* the x-slot, right? |
166 |
SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri) |
SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri) |
167 |
{ |
{ |
168 |
CHM_SP chxs = AS_CHM_SP__(x); |
CHM_SP chxs = AS_CHM_SP__(x); |
169 |
CHM_SP chxcp = cholmod_l_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c); |
CHM_SP chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c); |
170 |
int tr = asLogical(tri); |
int tr = asLogical(tri); |
171 |
R_CheckStack(); |
R_CheckStack(); |
172 |
|
|
176 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
177 |
} |
} |
178 |
|
|
179 |
|
// n.CMatrix --> [dli].CMatrix (not going through CHM!) |
180 |
|
SEXP nz_pattern_to_Csparse(SEXP x, SEXP res_kind) |
181 |
|
{ |
182 |
|
return nz2Csparse(x, asInteger(res_kind)); |
183 |
|
} |
184 |
|
// n.CMatrix --> [dli].CMatrix (not going through CHM!) |
185 |
|
SEXP nz2Csparse(SEXP x, enum x_slot_kind r_kind) |
186 |
|
{ |
187 |
|
const char *cl_x = class_P(x); |
188 |
|
if(cl_x[0] != 'n') error(_("not a 'n.CMatrix'")); |
189 |
|
if(cl_x[2] != 'C') error(_("not a CsparseMatrix")); |
190 |
|
int nnz = LENGTH(GET_SLOT(x, Matrix_iSym)); |
191 |
|
SEXP ans; |
192 |
|
char *ncl = strdup(cl_x); |
193 |
|
double *dx_x; int *ix_x; |
194 |
|
ncl[0] = (r_kind == x_double ? 'd' : |
195 |
|
(r_kind == x_logical ? 'l' : |
196 |
|
/* else (for now): r_kind == x_integer : */ 'i')); |
197 |
|
PROTECT(ans = NEW_OBJECT(MAKE_CLASS(ncl))); |
198 |
|
// create a correct 'x' slot: |
199 |
|
switch(r_kind) { |
200 |
|
int i; |
201 |
|
case x_double: // 'd' |
202 |
|
dx_x = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)); |
203 |
|
for (i=0; i < nnz; i++) dx_x[i] = 1.; |
204 |
|
break; |
205 |
|
case x_logical: // 'l' |
206 |
|
ix_x = LOGICAL(ALLOC_SLOT(ans, Matrix_xSym, LGLSXP, nnz)); |
207 |
|
for (i=0; i < nnz; i++) ix_x[i] = TRUE; |
208 |
|
break; |
209 |
|
case x_integer: // 'i' |
210 |
|
ix_x = INTEGER(ALLOC_SLOT(ans, Matrix_xSym, INTSXP, nnz)); |
211 |
|
for (i=0; i < nnz; i++) ix_x[i] = 1; |
212 |
|
break; |
213 |
|
|
214 |
|
default: |
215 |
|
error(_("nz2Csparse(): invalid/non-implemented r_kind = %d"), |
216 |
|
r_kind); |
217 |
|
} |
218 |
|
|
219 |
|
// now copy all other slots : |
220 |
|
slot_dup(ans, x, Matrix_iSym); |
221 |
|
slot_dup(ans, x, Matrix_pSym); |
222 |
|
slot_dup(ans, x, Matrix_DimSym); |
223 |
|
slot_dup(ans, x, Matrix_DimNamesSym); |
224 |
|
if(ncl[1] != 'g') { // symmetric or triangular ... |
225 |
|
slot_dup_if_has(ans, x, Matrix_uploSym); |
226 |
|
slot_dup_if_has(ans, x, Matrix_diagSym); |
227 |
|
} |
228 |
|
UNPROTECT(1); |
229 |
|
return ans; |
230 |
|
} |
231 |
|
|
232 |
SEXP Csparse_to_matrix(SEXP x) |
SEXP Csparse_to_matrix(SEXP x) |
233 |
{ |
{ |
234 |
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), |
235 |
1 /*do_free*/, GET_SLOT(x, Matrix_DimNamesSym)); |
1 /*do_free*/, GET_SLOT(x, Matrix_DimNamesSym)); |
236 |
} |
} |
237 |
|
|
238 |
SEXP Csparse_to_Tsparse(SEXP x, SEXP tri) |
SEXP Csparse_to_Tsparse(SEXP x, SEXP tri) |
239 |
{ |
{ |
240 |
CHM_SP chxs = AS_CHM_SP__(x); |
CHM_SP chxs = AS_CHM_SP__(x); |
241 |
CHM_TR chxt = cholmod_l_sparse_to_triplet(chxs, &c); |
CHM_TR chxt = cholmod_sparse_to_triplet(chxs, &c); |
242 |
int tr = asLogical(tri); |
int tr = asLogical(tri); |
243 |
int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
244 |
R_CheckStack(); |
R_CheckStack(); |
258 |
|
|
259 |
if (!(chx->stype)) |
if (!(chx->stype)) |
260 |
error(_("Nonsymmetric matrix in Csparse_symmetric_to_general")); |
error(_("Nonsymmetric matrix in Csparse_symmetric_to_general")); |
261 |
chgx = cholmod_l_copy(chx, /* stype: */ 0, chx->xtype, &c); |
chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c); |
262 |
/* xtype: pattern, "real", complex or .. */ |
/* xtype: pattern, "real", complex or .. */ |
263 |
return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
264 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
271 |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
272 |
R_CheckStack(); |
R_CheckStack(); |
273 |
|
|
274 |
chgx = cholmod_l_copy(chx, /* stype: */ uploT, chx->xtype, &c); |
chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c); |
275 |
/* xtype: pattern, "real", complex or .. */ |
/* xtype: pattern, "real", complex or .. */ |
276 |
return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
277 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
283 |
* since cholmod (& cs) lacks sparse 'int' matrices */ |
* since cholmod (& cs) lacks sparse 'int' matrices */ |
284 |
CHM_SP chx = AS_CHM_SP__(x); |
CHM_SP chx = AS_CHM_SP__(x); |
285 |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
286 |
CHM_SP chxt = cholmod_l_transpose(chx, chx->xtype, &c); |
CHM_SP chxt = cholmod_transpose(chx, chx->xtype, &c); |
287 |
SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp; |
SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp; |
288 |
int tr = asLogical(tri); |
int tr = asLogical(tri); |
289 |
R_CheckStack(); |
R_CheckStack(); |
302 |
CHM_SP |
CHM_SP |
303 |
cha = AS_CHM_SP(a), |
cha = AS_CHM_SP(a), |
304 |
chb = AS_CHM_SP(b), |
chb = AS_CHM_SP(b), |
305 |
chc = cholmod_l_ssmult(cha, chb, /*out_stype:*/ 0, |
chc = cholmod_ssmult(cha, chb, /*out_stype:*/ 0, |
306 |
/* values:= is_numeric (T/F) */ cha->xtype > 0, |
/* values:= is_numeric (T/F) */ cha->xtype > 0, |
307 |
/*out sorted:*/ 1, &c); |
/*out sorted:*/ 1, &c); |
308 |
const char *cl_a = class_P(a), *cl_b = class_P(b); |
const char *cl_a = class_P(a), *cl_b = class_P(b); |
353 |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
354 |
R_CheckStack(); |
R_CheckStack(); |
355 |
|
|
356 |
chTr = cholmod_l_transpose((tr) ? chb : cha, chb->xtype, &c); |
chTr = cholmod_transpose((tr) ? chb : cha, chb->xtype, &c); |
357 |
chc = cholmod_l_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb, |
chc = cholmod_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb, |
358 |
/*out_stype:*/ 0, cha->xtype, /*out sorted:*/ 1, &c); |
/*out_stype:*/ 0, cha->xtype, /*out sorted:*/ 1, &c); |
359 |
cholmod_l_free_sparse(&chTr, &c); |
cholmod_free_sparse(&chTr, &c); |
360 |
|
|
361 |
/* Preserve triangularity and unit-triangularity if appropriate; |
/* Preserve triangularity and unit-triangularity if appropriate; |
362 |
* see Csparse_Csparse_prod() for comments */ |
* see Csparse_Csparse_prod() for comments */ |
382 |
CHM_SP cha = AS_CHM_SP(a); |
CHM_SP cha = AS_CHM_SP(a); |
383 |
SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
384 |
CHM_DN chb = AS_CHM_DN(b_M); |
CHM_DN chb = AS_CHM_DN(b_M); |
385 |
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, |
386 |
chb->xtype, &c); |
chb->xtype, &c); |
387 |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
388 |
double one[] = {1,0}, zero[] = {0,0}; |
double one[] = {1,0}, zero[] = {0,0}; |
389 |
|
int nprot = 2; |
390 |
R_CheckStack(); |
R_CheckStack(); |
391 |
|
/* Tim Davis, please FIXME: currently (2010-11) *fails* when a is a pattern matrix:*/ |
392 |
cholmod_l_sdmult(cha, 0, one, zero, chb, chc, &c); |
if(cha->xtype == CHOLMOD_PATTERN) { |
393 |
|
/* warning(_("Csparse_dense_prod(): cholmod_sdmult() not yet implemented for pattern./ ngCMatrix" */ |
394 |
|
/* " --> slightly inefficient coercion")); */ |
395 |
|
|
396 |
|
// This *fails* to produce a CHOLMOD_REAL .. |
397 |
|
// CHM_SP chd = cholmod_l_copy(cha, cha->stype, CHOLMOD_REAL, &c); |
398 |
|
// --> use our Matrix-classes |
399 |
|
SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++; |
400 |
|
cha = AS_CHM_SP(da); |
401 |
|
} |
402 |
|
cholmod_sdmult(cha, 0, one, zero, chb, chc, &c); |
403 |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
404 |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); |
405 |
SET_VECTOR_ELT(dn, 1, |
SET_VECTOR_ELT(dn, 1, |
406 |
duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1))); |
duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1))); |
407 |
UNPROTECT(2); |
UNPROTECT(nprot); |
408 |
return chm_dense_to_SEXP(chc, 1, 0, dn); |
return chm_dense_to_SEXP(chc, 1, 0, dn); |
409 |
} |
} |
410 |
|
|
413 |
CHM_SP cha = AS_CHM_SP(a); |
CHM_SP cha = AS_CHM_SP(a); |
414 |
SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
415 |
CHM_DN chb = AS_CHM_DN(b_M); |
CHM_DN chb = AS_CHM_DN(b_M); |
416 |
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, |
417 |
chb->xtype, &c); |
chb->xtype, &c); |
418 |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); int nprot = 2; |
419 |
double one[] = {1,0}, zero[] = {0,0}; |
double one[] = {1,0}, zero[] = {0,0}; |
420 |
R_CheckStack(); |
R_CheckStack(); |
421 |
|
// -- see Csparse_dense_prod() above : |
422 |
cholmod_l_sdmult(cha, 1, one, zero, chb, chc, &c); |
if(cha->xtype == CHOLMOD_PATTERN) { |
423 |
|
SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++; |
424 |
|
cha = AS_CHM_SP(da); |
425 |
|
} |
426 |
|
cholmod_sdmult(cha, 1, one, zero, chb, chc, &c); |
427 |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
428 |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1))); |
duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1))); |
429 |
SET_VECTOR_ELT(dn, 1, |
SET_VECTOR_ELT(dn, 1, |
430 |
duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1))); |
duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1))); |
431 |
UNPROTECT(2); |
UNPROTECT(nprot); |
432 |
return chm_dense_to_SEXP(chc, 1, 0, dn); |
return chm_dense_to_SEXP(chc, 1, 0, dn); |
433 |
} |
} |
434 |
|
|
446 |
#endif |
#endif |
447 |
CHM_SP chcp, chxt, |
CHM_SP chcp, chxt, |
448 |
chx = (trip ? |
chx = (trip ? |
449 |
cholmod_l_triplet_to_sparse(cht, cht->nnz, &c) : |
cholmod_triplet_to_sparse(cht, cht->nnz, &c) : |
450 |
AS_CHM_SP(x)); |
AS_CHM_SP(x)); |
451 |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
452 |
R_CheckStack(); |
R_CheckStack(); |
453 |
|
|
454 |
if (!tr) chxt = cholmod_l_transpose(chx, chx->xtype, &c); |
if (!tr) chxt = cholmod_transpose(chx, chx->xtype, &c); |
455 |
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); |
456 |
if(!chcp) { |
if(!chcp) { |
457 |
UNPROTECT(1); |
UNPROTECT(1); |
458 |
error(_("Csparse_crossprod(): error return from cholmod_l_aat()")); |
error(_("Csparse_crossprod(): error return from cholmod_aat()")); |
459 |
} |
} |
460 |
cholmod_l_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c); |
cholmod_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c); |
461 |
chcp->stype = 1; |
chcp->stype = 1; |
462 |
if (trip) cholmod_l_free_sparse(&chx, &c); |
if (trip) cholmod_free_sparse(&chx, &c); |
463 |
if (!tr) cholmod_l_free_sparse(&chxt, &c); |
if (!tr) cholmod_free_sparse(&chxt, &c); |
464 |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
465 |
duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), |
duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), |
466 |
(tr) ? 0 : 1))); |
(tr) ? 0 : 1))); |
473 |
return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn); |
return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn); |
474 |
} |
} |
475 |
|
|
476 |
|
/* Csparse_drop(x, tol): drop entries with absolute value < tol, i.e, |
477 |
|
* at least all "explicit" zeros */ |
478 |
SEXP Csparse_drop(SEXP x, SEXP tol) |
SEXP Csparse_drop(SEXP x, SEXP tol) |
479 |
{ |
{ |
480 |
const char *cl = class_P(x); |
const char *cl = class_P(x); |
481 |
/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */ |
/* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */ |
482 |
int tr = (cl[1] == 't'); |
int tr = (cl[1] == 't'); |
483 |
CHM_SP chx = AS_CHM_SP__(x); |
CHM_SP chx = AS_CHM_SP__(x); |
484 |
CHM_SP ans = cholmod_l_copy(chx, chx->stype, chx->xtype, &c); |
CHM_SP ans = cholmod_copy(chx, chx->stype, chx->xtype, &c); |
485 |
double dtol = asReal(tol); |
double dtol = asReal(tol); |
486 |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
487 |
R_CheckStack(); |
R_CheckStack(); |
488 |
|
|
489 |
if(!cholmod_l_drop(dtol, ans, &c)) |
if(!cholmod_drop(dtol, ans, &c)) |
490 |
error(_("cholmod_l_drop() failed")); |
error(_("cholmod_drop() failed")); |
491 |
return chm_sparse_to_SEXP(ans, 1, |
return chm_sparse_to_SEXP(ans, 1, |
492 |
tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, |
tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, |
493 |
Rkind, tr ? diag_P(x) : "", |
Rkind, tr ? diag_P(x) : "", |
503 |
R_CheckStack(); |
R_CheckStack(); |
504 |
|
|
505 |
/* TODO: currently drops dimnames - and we fix at R level */ |
/* TODO: currently drops dimnames - and we fix at R level */ |
506 |
return chm_sparse_to_SEXP(cholmod_l_horzcat(chx, chy, 1, &c), |
return chm_sparse_to_SEXP(cholmod_horzcat(chx, chy, 1, &c), |
507 |
1, 0, Rkind, "", R_NilValue); |
1, 0, Rkind, "", R_NilValue); |
508 |
} |
} |
509 |
|
|
516 |
R_CheckStack(); |
R_CheckStack(); |
517 |
|
|
518 |
/* TODO: currently drops dimnames - and we fix at R level */ |
/* TODO: currently drops dimnames - and we fix at R level */ |
519 |
return chm_sparse_to_SEXP(cholmod_l_vertcat(chx, chy, 1, &c), |
return chm_sparse_to_SEXP(cholmod_vertcat(chx, chy, 1, &c), |
520 |
1, 0, Rkind, "", R_NilValue); |
1, 0, Rkind, "", R_NilValue); |
521 |
} |
} |
522 |
|
|
524 |
{ |
{ |
525 |
CHM_SP chx = AS_CHM_SP__(x); |
CHM_SP chx = AS_CHM_SP__(x); |
526 |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
527 |
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); |
528 |
R_CheckStack(); |
R_CheckStack(); |
529 |
|
|
530 |
return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", |
return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", |
542 |
} |
} |
543 |
else { /* unit triangular (diag='U'): "fill the diagonal" & diag:= "N" */ |
else { /* unit triangular (diag='U'): "fill the diagonal" & diag:= "N" */ |
544 |
CHM_SP chx = AS_CHM_SP__(x); |
CHM_SP chx = AS_CHM_SP__(x); |
545 |
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); |
546 |
double one[] = {1, 0}; |
double one[] = {1, 0}; |
547 |
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); |
548 |
int uploT = (*uplo_P(x) == 'U') ? 1 : -1; |
int uploT = (*uplo_P(x) == 'U') ? 1 : -1; |
549 |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
550 |
|
|
551 |
R_CheckStack(); |
R_CheckStack(); |
552 |
cholmod_l_free_sparse(&eye, &c); |
cholmod_free_sparse(&eye, &c); |
553 |
return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N", |
return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N", |
554 |
GET_SLOT(x, Matrix_DimNamesSym)); |
GET_SLOT(x, Matrix_DimNamesSym)); |
555 |
} |
} |
601 |
if (csize >= 0 && !isInteger(j)) |
if (csize >= 0 && !isInteger(j)) |
602 |
error(_("Index j must be NULL or integer")); |
error(_("Index j must be NULL or integer")); |
603 |
|
|
604 |
if (chx->stype) /* symmetricMatrix */ |
if (!chx->stype) {/* non-symmetric Matrix */ |
605 |
/* for now, cholmod_submatrix() only accepts "generalMatrix" */ |
return chm_sparse_to_SEXP(cholmod_submatrix(chx, |
|
chx = cholmod_l_copy(chx, /* stype: */ 0, chx->xtype, &c); |
|
|
|
|
|
return chm_sparse_to_SEXP(cholmod_l_submatrix(chx, |
|
606 |
(rsize < 0) ? NULL : INTEGER(i), rsize, |
(rsize < 0) ? NULL : INTEGER(i), rsize, |
607 |
(csize < 0) ? NULL : INTEGER(j), csize, |
(csize < 0) ? NULL : INTEGER(j), csize, |
608 |
TRUE, TRUE, &c), |
TRUE, TRUE, &c), |
609 |
1, 0, Rkind, "", |
1, 0, Rkind, "", |
610 |
/* FIXME: drops dimnames */ R_NilValue); |
/* FIXME: drops dimnames */ R_NilValue); |
611 |
} |
} |
612 |
|
/* for now, cholmod_submatrix() only accepts "generalMatrix" */ |
613 |
|
CHM_SP tmp = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c); |
614 |
|
CHM_SP ans = cholmod_submatrix(tmp, |
615 |
|
(rsize < 0) ? NULL : INTEGER(i), rsize, |
616 |
|
(csize < 0) ? NULL : INTEGER(j), csize, |
617 |
|
TRUE, TRUE, &c); |
618 |
|
cholmod_free_sparse(&tmp, &c); |
619 |
|
return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue); |
620 |
|
} |
621 |
|
|
622 |
|
#define _d_Csp_ |
623 |
|
#include "t_Csparse_subassign.c" |
624 |
|
|
625 |
|
#define _l_Csp_ |
626 |
|
#include "t_Csparse_subassign.c" |
627 |
|
|
628 |
|
#define _i_Csp_ |
629 |
|
#include "t_Csparse_subassign.c" |
630 |
|
|
631 |
|
#define _n_Csp_ |
632 |
|
#include "t_Csparse_subassign.c" |
633 |
|
|
634 |
|
#define _z_Csp_ |
635 |
|
#include "t_Csparse_subassign.c" |
636 |
|
|
637 |
|
|
638 |
|
|
639 |
SEXP Csparse_MatrixMarket(SEXP x, SEXP fname) |
SEXP Csparse_MatrixMarket(SEXP x, SEXP fname) |
640 |
{ |
{ |
643 |
if (!f) |
if (!f) |
644 |
error(_("failure to open file \"%s\" for writing"), |
error(_("failure to open file \"%s\" for writing"), |
645 |
CHAR(asChar(fname))); |
CHAR(asChar(fname))); |
646 |
if (!cholmod_l_write_sparse(f, AS_CHM_SP(x), |
if (!cholmod_write_sparse(f, AS_CHM_SP(x), |
647 |
(CHM_SP)NULL, (char*) NULL, &c)) |
(CHM_SP)NULL, (char*) NULL, &c)) |
648 |
error(_("cholmod_l_write_sparse returned error code")); |
error(_("cholmod_write_sparse returned error code")); |
649 |
fclose(f); |
fclose(f); |
650 |
return R_NilValue; |
return R_NilValue; |
651 |
} |
} |
859 |
if (cls[1] != 'g') |
if (cls[1] != 'g') |
860 |
error(_("Only 'g'eneral sparse matrix types allowed")); |
error(_("Only 'g'eneral sparse matrix types allowed")); |
861 |
/* allocate and populate the triplet */ |
/* allocate and populate the triplet */ |
862 |
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, |
863 |
xtype, &c); |
xtype, &c); |
864 |
T->x = x; |
T->x = x; |
865 |
tri = (int*)T->i; |
tri = (int*)T->i; |
869 |
trj[ii] = j[ii] - ((!mj && index1) ? 1 : 0); |
trj[ii] = j[ii] - ((!mj && index1) ? 1 : 0); |
870 |
} |
} |
871 |
/* create the cholmod_sparse structure */ |
/* create the cholmod_sparse structure */ |
872 |
A = cholmod_l_triplet_to_sparse(T, nnz, &c); |
A = cholmod_triplet_to_sparse(T, nnz, &c); |
873 |
cholmod_l_free_triplet(&T, &c); |
cholmod_free_triplet(&T, &c); |
874 |
/* copy the information to the SEXP */ |
/* copy the information to the SEXP */ |
875 |
ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cls))); |
ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cls))); |
876 |
/* 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 */ |
877 |
/* allocate and copy common slots */ |
/* allocate and copy common slots */ |
878 |
nnz = cholmod_l_nnz(A, &c); |
nnz = cholmod_nnz(A, &c); |
879 |
dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)); |
dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)); |
880 |
dims[0] = A->nrow; dims[1] = A->ncol; |
dims[0] = A->nrow; dims[1] = A->ncol; |
881 |
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); |
887 |
case 'l': |
case 'l': |
888 |
error(_("code not yet written for cls = \"lgCMatrix\"")); |
error(_("code not yet written for cls = \"lgCMatrix\"")); |
889 |
} |
} |
890 |
cholmod_l_free_sparse(&A, &c); |
/* FIXME: dimnames are *NOT* put there yet (if non-NULL) */ |
891 |
|
cholmod_free_sparse(&A, &c); |
892 |
UNPROTECT(1); |
UNPROTECT(1); |
893 |
return ans; |
return ans; |
894 |
} |
} |