SCM Repository
Annotation of /pkg/src/Csparse.c
Parent Directory
| 
Revision Log
Revision 2299 - (view) (download) (as text)
| 1 : | bates | 1218 | /* Sparse matrices in compressed column-oriented form */ |
| 2 : | bates | 922 | #include "Csparse.h" |
| 3 : | maechler | 2120 | #include "Tsparse.h" |
| 4 : | bates | 922 | #include "chm_common.h" |
| 5 : | |||
| 6 : | mmaechler | 2279 | /** "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 */ | ||
| 10 : | SEXP pslot = GET_SLOT(x, Matrix_pSym), | ||
| 11 : | islot = GET_SLOT(x, Matrix_iSym); | ||
| 12 : | 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 : | bates | 922 | SEXP Csparse_validate(SEXP x) |
| 36 : | { | ||
| 37 : | maechler | 1575 | /* NB: we do *NOT* check a potential 'x' slot here, at all */ |
| 38 : | bates | 922 | SEXP pslot = GET_SLOT(x, Matrix_pSym), |
| 39 : | islot = GET_SLOT(x, Matrix_iSym); | ||
| 40 : | maechler | 1893 | Rboolean sorted, strictly; |
| 41 : | int j, k, | ||
| 42 : | bates | 922 | *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), |
| 43 : | maechler | 1660 | nrow = dims[0], |
| 44 : | ncol = dims[1], | ||
| 45 : | maechler | 1654 | *xp = INTEGER(pslot), |
| 46 : | bates | 922 | *xi = INTEGER(islot); |
| 47 : | |||
| 48 : | maechler | 1654 | if (length(pslot) != dims[1] + 1) |
| 49 : | return mkString(_("slot p must have length = ncol(.) + 1")); | ||
| 50 : | bates | 922 | if (xp[0] != 0) |
| 51 : | return mkString(_("first element of slot p must be zero")); | ||
| 52 : | maechler | 1893 | if (length(islot) < xp[ncol]) /* allow larger slots from over-allocation!*/ |
| 53 : | bates | 1555 | return |
| 54 : | mkString(_("last element of slot p must match length of slots i and x")); | ||
| 55 : | mmaechler | 2236 | for (j = 0; j < xp[ncol]; j++) { |
| 56 : | bates | 1555 | if (xi[j] < 0 || xi[j] >= nrow) |
| 57 : | return mkString(_("all row indices must be between 0 and nrow-1")); | ||
| 58 : | } | ||
| 59 : | maechler | 1893 | sorted = TRUE; strictly = TRUE; |
| 60 : | bates | 922 | for (j = 0; j < ncol; j++) { |
| 61 : | mmaechler | 2236 | if (xp[j] > xp[j + 1]) |
| 62 : | bates | 922 | return mkString(_("slot p must be non-decreasing")); |
| 63 : | mmaechler | 2236 | if(sorted) /* only act if >= 2 entries in column j : */ |
| 64 : | maechler | 1893 | for (k = xp[j] + 1; k < xp[j + 1]; k++) { |
| 65 : | if (xi[k] < xi[k - 1]) | ||
| 66 : | sorted = FALSE; | ||
| 67 : | else if (xi[k] == xi[k - 1]) | ||
| 68 : | strictly = FALSE; | ||
| 69 : | } | ||
| 70 : | bates | 922 | } |
| 71 : | maechler | 1654 | if (!sorted) { |
| 72 : | mmaechler | 2235 | CHM_SP chx = (CHM_SP) alloca(sizeof(cholmod_sparse)); |
| 73 : | maechler | 1960 | R_CheckStack(); |
| 74 : | dmbates | 2299 | as_cholmod_sparse(chx, x, FALSE, TRUE); /* includes cholmod_l_sort() ! */ |
| 75 : | mmaechler | 2235 | /* as chx = AS_CHM_SP__(x) but ^^^^ sorting x in_place (no copying)*/ |
| 76 : | maechler | 1960 | |
| 77 : | maechler | 1893 | /* Now re-check that row indices are *strictly* increasing |
| 78 : | * (and not just increasing) within each column : */ | ||
| 79 : | for (j = 0; j < ncol; j++) { | ||
| 80 : | for (k = xp[j] + 1; k < xp[j + 1]; k++) | ||
| 81 : | if (xi[k] == xi[k - 1]) | ||
| 82 : | dmbates | 2299 | return mkString(_("slot i is not *strictly* increasing inside a column (even after cholmod_l_sort)")); |
| 83 : | maechler | 1893 | } |
| 84 : | |||
| 85 : | } else if(!strictly) { /* sorted, but not strictly */ | ||
| 86 : | return mkString(_("slot i is not *strictly* increasing inside a column")); | ||
| 87 : | maechler | 1654 | } |
| 88 : | bates | 922 | return ScalarLogical(1); |
| 89 : | } | ||
| 90 : | |||
| 91 : | maechler | 1968 | SEXP Rsparse_validate(SEXP x) |
| 92 : | { | ||
| 93 : | /* NB: we do *NOT* check a potential 'x' slot here, at all */ | ||
| 94 : | SEXP pslot = GET_SLOT(x, Matrix_pSym), | ||
| 95 : | jslot = GET_SLOT(x, Matrix_jSym); | ||
| 96 : | Rboolean sorted, strictly; | ||
| 97 : | int i, k, | ||
| 98 : | *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), | ||
| 99 : | nrow = dims[0], | ||
| 100 : | ncol = dims[1], | ||
| 101 : | *xp = INTEGER(pslot), | ||
| 102 : | *xj = INTEGER(jslot); | ||
| 103 : | |||
| 104 : | if (length(pslot) != dims[0] + 1) | ||
| 105 : | return mkString(_("slot p must have length = nrow(.) + 1")); | ||
| 106 : | if (xp[0] != 0) | ||
| 107 : | return mkString(_("first element of slot p must be zero")); | ||
| 108 : | if (length(jslot) < xp[nrow]) /* allow larger slots from over-allocation!*/ | ||
| 109 : | return | ||
| 110 : | mkString(_("last element of slot p must match length of slots j and x")); | ||
| 111 : | for (i = 0; i < length(jslot); i++) { | ||
| 112 : | if (xj[i] < 0 || xj[i] >= ncol) | ||
| 113 : | return mkString(_("all column indices must be between 0 and ncol-1")); | ||
| 114 : | } | ||
| 115 : | sorted = TRUE; strictly = TRUE; | ||
| 116 : | for (i = 0; i < nrow; i++) { | ||
| 117 : | if (xp[i] > xp[i+1]) | ||
| 118 : | return mkString(_("slot p must be non-decreasing")); | ||
| 119 : | if(sorted) | ||
| 120 : | for (k = xp[i] + 1; k < xp[i + 1]; k++) { | ||
| 121 : | if (xj[k] < xj[k - 1]) | ||
| 122 : | sorted = FALSE; | ||
| 123 : | else if (xj[k] == xj[k - 1]) | ||
| 124 : | strictly = FALSE; | ||
| 125 : | } | ||
| 126 : | } | ||
| 127 : | if (!sorted) | ||
| 128 : | dmbates | 2299 | /* cannot easily use cholmod_l_sort(.) ... -> "error out" :*/ |
| 129 : | maechler | 1968 | return mkString(_("slot j is not increasing inside a column")); |
| 130 : | else if(!strictly) /* sorted, but not strictly */ | ||
| 131 : | return mkString(_("slot j is not *strictly* increasing inside a column")); | ||
| 132 : | |||
| 133 : | return ScalarLogical(1); | ||
| 134 : | } | ||
| 135 : | |||
| 136 : | |||
| 137 : | maechler | 1751 | /* Called from ../R/Csparse.R : */ |
| 138 : | /* Can only return [dln]geMatrix (no symm/triang); | ||
| 139 : | * FIXME: replace by non-CHOLMOD code ! */ | ||
| 140 : | bates | 1059 | SEXP Csparse_to_dense(SEXP x) |
| 141 : | { | ||
| 142 : | mmaechler | 2223 | CHM_SP chxs = AS_CHM_SP__(x); |
| 143 : | maechler | 1751 | /* This loses the symmetry property, since cholmod_dense has none, |
| 144 : | * BUT, much worse (FIXME!), it also transforms CHOLMOD_PATTERN ("n") matrices | ||
| 145 : | * to numeric (CHOLMOD_REAL) ones : */ | ||
| 146 : | dmbates | 2299 | CHM_DN chxd = cholmod_l_sparse_to_dense(chxs, &c); |
| 147 : | maechler | 1751 | int Rkind = (chxs->xtype == CHOLMOD_PATTERN)? -1 : Real_kind(x); |
| 148 : | maechler | 1960 | R_CheckStack(); |
| 149 : | bates | 1059 | |
| 150 : | maechler | 1736 | return chm_dense_to_SEXP(chxd, 1, Rkind, GET_SLOT(x, Matrix_DimNamesSym)); |
| 151 : | bates | 1059 | } |
| 152 : | |||
| 153 : | maechler | 1548 | SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri) |
| 154 : | bates | 1371 | { |
| 155 : | mmaechler | 2223 | CHM_SP chxs = AS_CHM_SP__(x); |
| 156 : | dmbates | 2299 | CHM_SP chxcp = cholmod_l_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c); |
| 157 : | bates | 1867 | int tr = asLogical(tri); |
| 158 : | maechler | 1960 | R_CheckStack(); |
| 159 : | bates | 1371 | |
| 160 : | maechler | 1960 | return chm_sparse_to_SEXP(chxcp, 1/*do_free*/, |
| 161 : | bates | 1867 | tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, |
| 162 : | 0, tr ? diag_P(x) : "", | ||
| 163 : | maechler | 1548 | GET_SLOT(x, Matrix_DimNamesSym)); |
| 164 : | bates | 1371 | } |
| 165 : | |||
| 166 : | bates | 1366 | SEXP Csparse_to_matrix(SEXP x) |
| 167 : | bates | 922 | { |
| 168 : | dmbates | 2299 | return chm_dense_to_matrix(cholmod_l_sparse_to_dense(AS_CHM_SP__(x), &c), |
| 169 : | maechler | 1960 | 1 /*do_free*/, GET_SLOT(x, Matrix_DimNamesSym)); |
| 170 : | bates | 1366 | } |
| 171 : | |||
| 172 : | SEXP Csparse_to_Tsparse(SEXP x, SEXP tri) | ||
| 173 : | { | ||
| 174 : | mmaechler | 2223 | CHM_SP chxs = AS_CHM_SP__(x); |
| 175 : | dmbates | 2299 | CHM_TR chxt = cholmod_l_sparse_to_triplet(chxs, &c); |
| 176 : | bates | 1867 | int tr = asLogical(tri); |
| 177 : | maechler | 1736 | int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
| 178 : | maechler | 1960 | R_CheckStack(); |
| 179 : | bates | 922 | |
| 180 : | bates | 1867 | return chm_triplet_to_SEXP(chxt, 1, |
| 181 : | tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, | ||
| 182 : | Rkind, tr ? diag_P(x) : "", | ||
| 183 : | bates | 1371 | GET_SLOT(x, Matrix_DimNamesSym)); |
| 184 : | bates | 922 | } |
| 185 : | |||
| 186 : | bates | 1448 | /* this used to be called sCMatrix_to_gCMatrix(..) [in ./dsCMatrix.c ]: */ |
| 187 : | bates | 1371 | SEXP Csparse_symmetric_to_general(SEXP x) |
| 188 : | { | ||
| 189 : | mmaechler | 2223 | CHM_SP chx = AS_CHM_SP__(x), chgx; |
| 190 : | maechler | 1736 | int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
| 191 : | maechler | 1960 | R_CheckStack(); |
| 192 : | bates | 1371 | |
| 193 : | if (!(chx->stype)) | ||
| 194 : | maechler | 1548 | error(_("Nonsymmetric matrix in Csparse_symmetric_to_general")); |
| 195 : | dmbates | 2299 | chgx = cholmod_l_copy(chx, /* stype: */ 0, chx->xtype, &c); |
| 196 : | maechler | 1375 | /* xtype: pattern, "real", complex or .. */ |
| 197 : | maechler | 1548 | return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", |
| 198 : | bates | 1371 | GET_SLOT(x, Matrix_DimNamesSym)); |
| 199 : | } | ||
| 200 : | |||
| 201 : | maechler | 1618 | SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo) |
| 202 : | maechler | 1598 | { |
| 203 : | mmaechler | 2223 | CHM_SP chx = AS_CHM_SP__(x), chgx; |
| 204 : | maechler | 2113 | int uploT = (*CHAR(STRING_ELT(uplo,0)) == 'U') ? 1 : -1; |
| 205 : | maechler | 1736 | int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
| 206 : | maechler | 1960 | R_CheckStack(); |
| 207 : | maechler | 1598 | |
| 208 : | dmbates | 2299 | chgx = cholmod_l_copy(chx, /* stype: */ uploT, chx->xtype, &c); |
| 209 : | maechler | 1598 | /* xtype: pattern, "real", complex or .. */ |
| 210 : | return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", | ||
| 211 : | GET_SLOT(x, Matrix_DimNamesSym)); | ||
| 212 : | } | ||
| 213 : | |||
| 214 : | bates | 1369 | SEXP Csparse_transpose(SEXP x, SEXP tri) |
| 215 : | bates | 922 | { |
| 216 : | maechler | 1921 | /* TODO: lgCMatrix & igC* currently go via double prec. cholmod - |
| 217 : | * since cholmod (& cs) lacks sparse 'int' matrices */ | ||
| 218 : | mmaechler | 2223 | CHM_SP chx = AS_CHM_SP__(x); |
| 219 : | maechler | 1736 | int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
| 220 : | dmbates | 2299 | CHM_SP chxt = cholmod_l_transpose(chx, chx->xtype, &c); |
| 221 : | bates | 1366 | SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp; |
| 222 : | bates | 1867 | int tr = asLogical(tri); |
| 223 : | maechler | 1960 | R_CheckStack(); |
| 224 : | bates | 1369 | |
| 225 : | bates | 1366 | tmp = VECTOR_ELT(dn, 0); /* swap the dimnames */ |
| 226 : | SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1)); | ||
| 227 : | SET_VECTOR_ELT(dn, 1, tmp); | ||
| 228 : | UNPROTECT(1); | ||
| 229 : | bates | 1867 | return chm_sparse_to_SEXP(chxt, 1, /* SWAP 'uplo' for triangular */ |
| 230 : | tr ? ((*uplo_P(x) == 'U') ? -1 : 1) : 0, | ||
| 231 : | Rkind, tr ? diag_P(x) : "", dn); | ||
| 232 : | bates | 922 | } |
| 233 : | |||
| 234 : | SEXP Csparse_Csparse_prod(SEXP a, SEXP b) | ||
| 235 : | { | ||
| 236 : | maechler | 2120 | CHM_SP |
| 237 : | mmaechler | 2223 | cha = AS_CHM_SP(a), |
| 238 : | chb = AS_CHM_SP(b), | ||
| 239 : | dmbates | 2299 | chc = cholmod_l_ssmult(cha, chb, /*out_stype:*/ 0, |
| 240 : | maechler | 2125 | cha->xtype, /*out sorted:*/ 1, &c); |
| 241 : | const char *cl_a = class_P(a), *cl_b = class_P(b); | ||
| 242 : | char diag[] = {'\0', '\0'}; | ||
| 243 : | int uploT = 0; | ||
| 244 : | bates | 1366 | SEXP dn = allocVector(VECSXP, 2); |
| 245 : | maechler | 1960 | R_CheckStack(); |
| 246 : | bates | 922 | |
| 247 : | maechler | 2125 | /* Preserve triangularity and even unit-triangularity if appropriate. |
| 248 : | * Note that in that case, the multiplication itself should happen | ||
| 249 : | * faster. But there's no support for that in CHOLMOD */ | ||
| 250 : | |||
| 251 : | /* UGLY hack -- rather should have (fast!) C-level version of | ||
| 252 : | * is(a, "triangularMatrix") etc */ | ||
| 253 : | if (cl_a[1] == 't' && cl_b[1] == 't') | ||
| 254 : | /* FIXME: fails for "Cholesky","BunchKaufmann"..*/ | ||
| 255 : | if(*uplo_P(a) == *uplo_P(b)) { /* both upper, or both lower tri. */ | ||
| 256 : | uploT = (*uplo_P(a) == 'U') ? 1 : -1; | ||
| 257 : | if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */ | ||
| 258 : | /* "remove the diagonal entries": */ | ||
| 259 : | chm_diagN2U(chc, uploT, /* do_realloc */ FALSE); | ||
| 260 : | diag[0]= 'U'; | ||
| 261 : | } | ||
| 262 : | else diag[0]= 'N'; | ||
| 263 : | } | ||
| 264 : | bates | 1366 | SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
| 265 : | duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); | ||
| 266 : | SET_VECTOR_ELT(dn, 1, | ||
| 267 : | duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1))); | ||
| 268 : | maechler | 2125 | return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn); |
| 269 : | bates | 922 | } |
| 270 : | |||
| 271 : | maechler | 1659 | SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b, SEXP trans) |
| 272 : | bates | 1657 | { |
| 273 : | maechler | 1659 | int tr = asLogical(trans); |
| 274 : | maechler | 2120 | CHM_SP |
| 275 : | mmaechler | 2223 | cha = AS_CHM_SP(a), |
| 276 : | chb = AS_CHM_SP(b), | ||
| 277 : | maechler | 2120 | chTr, chc; |
| 278 : | maechler | 2125 | const char *cl_a = class_P(a), *cl_b = class_P(b); |
| 279 : | char diag[] = {'\0', '\0'}; | ||
| 280 : | int uploT = 0; | ||
| 281 : | bates | 1657 | SEXP dn = allocVector(VECSXP, 2); |
| 282 : | maechler | 1960 | R_CheckStack(); |
| 283 : | bates | 1657 | |
| 284 : | dmbates | 2299 | chTr = cholmod_l_transpose((tr) ? chb : cha, chb->xtype, &c); |
| 285 : | chc = cholmod_l_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb, | ||
| 286 : | maechler | 2125 | /*out_stype:*/ 0, cha->xtype, /*out sorted:*/ 1, &c); |
| 287 : | dmbates | 2299 | cholmod_l_free_sparse(&chTr, &c); |
| 288 : | maechler | 1659 | |
| 289 : | maechler | 2125 | /* Preserve triangularity and unit-triangularity if appropriate; |
| 290 : | * see Csparse_Csparse_prod() for comments */ | ||
| 291 : | if (cl_a[1] == 't' && cl_b[1] == 't') | ||
| 292 : | if(*uplo_P(a) != *uplo_P(b)) { /* one 'U', the other 'L' */ | ||
| 293 : | uploT = (*uplo_P(b) == 'U') ? 1 : -1; | ||
| 294 : | if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */ | ||
| 295 : | chm_diagN2U(chc, uploT, /* do_realloc */ FALSE); | ||
| 296 : | diag[0]= 'U'; | ||
| 297 : | } | ||
| 298 : | else diag[0]= 'N'; | ||
| 299 : | } | ||
| 300 : | |||
| 301 : | bates | 1657 | SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
| 302 : | maechler | 1659 | duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), (tr) ? 0 : 1))); |
| 303 : | bates | 1657 | SET_VECTOR_ELT(dn, 1, |
| 304 : | maechler | 1659 | duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), (tr) ? 0 : 1))); |
| 305 : | maechler | 2125 | return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn); |
| 306 : | bates | 1657 | } |
| 307 : | |||
| 308 : | bates | 922 | SEXP Csparse_dense_prod(SEXP a, SEXP b) |
| 309 : | { | ||
| 310 : | mmaechler | 2223 | CHM_SP cha = AS_CHM_SP(a); |
| 311 : | maechler | 1660 | SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
| 312 : | maechler | 1960 | CHM_DN chb = AS_CHM_DN(b_M); |
| 313 : | dmbates | 2299 | CHM_DN chc = cholmod_l_allocate_dense(cha->nrow, chb->ncol, cha->nrow, |
| 314 : | maechler | 1960 | chb->xtype, &c); |
| 315 : | SEXP dn = PROTECT(allocVector(VECSXP, 2)); | ||
| 316 : | double one[] = {1,0}, zero[] = {0,0}; | ||
| 317 : | R_CheckStack(); | ||
| 318 : | bates | 922 | |
| 319 : | dmbates | 2299 | cholmod_l_sdmult(cha, 0, one, zero, chb, chc, &c); |
| 320 : | maechler | 1660 | SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
| 321 : | duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); | ||
| 322 : | SET_VECTOR_ELT(dn, 1, | ||
| 323 : | duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1))); | ||
| 324 : | maechler | 1960 | UNPROTECT(2); |
| 325 : | maechler | 1660 | return chm_dense_to_SEXP(chc, 1, 0, dn); |
| 326 : | bates | 922 | } |
| 327 : | maechler | 925 | |
| 328 : | bates | 1067 | SEXP Csparse_dense_crossprod(SEXP a, SEXP b) |
| 329 : | { | ||
| 330 : | mmaechler | 2223 | CHM_SP cha = AS_CHM_SP(a); |
| 331 : | maechler | 1660 | SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b)); |
| 332 : | maechler | 1960 | CHM_DN chb = AS_CHM_DN(b_M); |
| 333 : | dmbates | 2299 | CHM_DN chc = cholmod_l_allocate_dense(cha->ncol, chb->ncol, cha->ncol, |
| 334 : | maechler | 1960 | chb->xtype, &c); |
| 335 : | SEXP dn = PROTECT(allocVector(VECSXP, 2)); | ||
| 336 : | double one[] = {1,0}, zero[] = {0,0}; | ||
| 337 : | R_CheckStack(); | ||
| 338 : | bates | 1067 | |
| 339 : | dmbates | 2299 | cholmod_l_sdmult(cha, 1, one, zero, chb, chc, &c); |
| 340 : | maechler | 1660 | SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
| 341 : | duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1))); | ||
| 342 : | SET_VECTOR_ELT(dn, 1, | ||
| 343 : | duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1))); | ||
| 344 : | maechler | 1960 | UNPROTECT(2); |
| 345 : | maechler | 1660 | return chm_dense_to_SEXP(chc, 1, 0, dn); |
| 346 : | bates | 1067 | } |
| 347 : | |||
| 348 : | maechler | 2125 | /* Computes x'x or x x' -- *also* for Tsparse (triplet = TRUE) |
| 349 : | see Csparse_Csparse_crossprod above for x'y and x y' */ | ||
| 350 : | bates | 928 | SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet) |
| 351 : | bates | 922 | { |
| 352 : | maechler | 957 | int trip = asLogical(triplet), |
| 353 : | tr = asLogical(trans); /* gets reversed because _aat is tcrossprod */ | ||
| 354 : | mmaechler | 2223 | CHM_TR cht = trip ? AS_CHM_TR(x) : (CHM_TR) NULL; |
| 355 : | maechler | 1960 | CHM_SP chcp, chxt, |
| 356 : | maechler | 2120 | chx = (trip ? |
| 357 : | dmbates | 2299 | cholmod_l_triplet_to_sparse(cht, cht->nnz, &c) : |
| 358 : | mmaechler | 2223 | AS_CHM_SP(x)); |
| 359 : | bates | 1366 | SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
| 360 : | maechler | 1960 | R_CheckStack(); |
| 361 : | bates | 922 | |
| 362 : | dmbates | 2299 | if (!tr) chxt = cholmod_l_transpose(chx, chx->xtype, &c); |
| 363 : | chcp = cholmod_l_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c); | ||
| 364 : | maechler | 2120 | if(!chcp) { |
| 365 : | UNPROTECT(1); | ||
| 366 : | dmbates | 2299 | error(_("Csparse_crossprod(): error return from cholmod_l_aat()")); |
| 367 : | maechler | 2120 | } |
| 368 : | dmbates | 2299 | cholmod_l_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c); |
| 369 : | bates | 1360 | chcp->stype = 1; |
| 370 : | dmbates | 2299 | if (trip) cholmod_l_free_sparse(&chx, &c); |
| 371 : | if (!tr) cholmod_l_free_sparse(&chxt, &c); | ||
| 372 : | maechler | 1960 | SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
| 373 : | bates | 1366 | duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), |
| 374 : | maechler | 1660 | (tr) ? 0 : 1))); |
| 375 : | bates | 1366 | SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(dn, 0))); |
| 376 : | UNPROTECT(1); | ||
| 377 : | maechler | 1548 | return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn); |
| 378 : | bates | 922 | } |
| 379 : | bates | 923 | |
| 380 : | maechler | 1618 | SEXP Csparse_drop(SEXP x, SEXP tol) |
| 381 : | { | ||
| 382 : | mmaechler | 2175 | const char *cl = class_P(x); |
| 383 : | /* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */ | ||
| 384 : | int tr = (cl[1] == 't'); | ||
| 385 : | mmaechler | 2223 | CHM_SP chx = AS_CHM_SP__(x); |
| 386 : | dmbates | 2299 | CHM_SP ans = cholmod_l_copy(chx, chx->stype, chx->xtype, &c); |
| 387 : | maechler | 1618 | double dtol = asReal(tol); |
| 388 : | maechler | 1736 | int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
| 389 : | maechler | 1960 | R_CheckStack(); |
| 390 : | maechler | 1618 | |
| 391 : | dmbates | 2299 | if(!cholmod_l_drop(dtol, ans, &c)) |
| 392 : | error(_("cholmod_l_drop() failed")); | ||
| 393 : | mmaechler | 2175 | return chm_sparse_to_SEXP(ans, 1, |
| 394 : | tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, | ||
| 395 : | Rkind, tr ? diag_P(x) : "", | ||
| 396 : | maechler | 1736 | GET_SLOT(x, Matrix_DimNamesSym)); |
| 397 : | maechler | 1618 | } |
| 398 : | |||
| 399 : | bates | 1218 | SEXP Csparse_horzcat(SEXP x, SEXP y) |
| 400 : | { | ||
| 401 : | mmaechler | 2223 | CHM_SP chx = AS_CHM_SP__(x), chy = AS_CHM_SP__(y); |
| 402 : | maechler | 1548 | int Rkind = 0; /* only for "d" - FIXME */ |
| 403 : | maechler | 1960 | R_CheckStack(); |
| 404 : | maechler | 1375 | |
| 405 : | bates | 1366 | /* FIXME: currently drops dimnames */ |
| 406 : | dmbates | 2299 | return chm_sparse_to_SEXP(cholmod_l_horzcat(chx, chy, 1, &c), |
| 407 : | maechler | 1960 | 1, 0, Rkind, "", R_NilValue); |
| 408 : | bates | 1218 | } |
| 409 : | |||
| 410 : | SEXP Csparse_vertcat(SEXP x, SEXP y) | ||
| 411 : | { | ||
| 412 : | mmaechler | 2223 | CHM_SP chx = AS_CHM_SP__(x), chy = AS_CHM_SP__(y); |
| 413 : | maechler | 1548 | int Rkind = 0; /* only for "d" - FIXME */ |
| 414 : | maechler | 1960 | R_CheckStack(); |
| 415 : | maechler | 1375 | |
| 416 : | bates | 1366 | /* FIXME: currently drops dimnames */ |
| 417 : | dmbates | 2299 | return chm_sparse_to_SEXP(cholmod_l_vertcat(chx, chy, 1, &c), |
| 418 : | maechler | 1960 | 1, 0, Rkind, "", R_NilValue); |
| 419 : | bates | 1218 | } |
| 420 : | bates | 1265 | |
| 421 : | SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2) | ||
| 422 : | { | ||
| 423 : | mmaechler | 2223 | CHM_SP chx = AS_CHM_SP__(x); |
| 424 : | maechler | 1736 | int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
| 425 : | dmbates | 2299 | CHM_SP ans = cholmod_l_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c); |
| 426 : | maechler | 1960 | R_CheckStack(); |
| 427 : | bates | 1265 | |
| 428 : | maechler | 1736 | return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", |
| 429 : | GET_SLOT(x, Matrix_DimNamesSym)); | ||
| 430 : | bates | 1265 | } |
| 431 : | bates | 1366 | |
| 432 : | SEXP Csparse_diagU2N(SEXP x) | ||
| 433 : | { | ||
| 434 : | maechler | 2120 | const char *cl = class_P(x); |
| 435 : | /* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */ | ||
| 436 : | if (cl[1] != 't' || *diag_P(x) != 'U') { | ||
| 437 : | maechler | 2125 | /* "trivially fast" when not triangular (<==> no 'diag' slot), |
| 438 : | or not *unit* triangular */ | ||
| 439 : | maechler | 1708 | return (x); |
| 440 : | } | ||
| 441 : | maechler | 2125 | else { /* unit triangular (diag='U'): "fill the diagonal" & diag:= "N" */ |
| 442 : | mmaechler | 2223 | CHM_SP chx = AS_CHM_SP__(x); |
| 443 : | dmbates | 2299 | CHM_SP eye = cholmod_l_speye(chx->nrow, chx->ncol, chx->xtype, &c); |
| 444 : | maechler | 1708 | double one[] = {1, 0}; |
| 445 : | dmbates | 2299 | CHM_SP ans = cholmod_l_add(chx, eye, one, one, TRUE, TRUE, &c); |
| 446 : | maechler | 1710 | int uploT = (*uplo_P(x) == 'U') ? 1 : -1; |
| 447 : | maechler | 1736 | int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
| 448 : | bates | 1366 | |
| 449 : | maechler | 1960 | R_CheckStack(); |
| 450 : | dmbates | 2299 | cholmod_l_free_sparse(&eye, &c); |
| 451 : | maechler | 1708 | return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N", |
| 452 : | maechler | 1736 | GET_SLOT(x, Matrix_DimNamesSym)); |
| 453 : | maechler | 1708 | } |
| 454 : | bates | 1366 | } |
| 455 : | |||
| 456 : | maechler | 2125 | SEXP Csparse_diagN2U(SEXP x) |
| 457 : | { | ||
| 458 : | const char *cl = class_P(x); | ||
| 459 : | /* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */ | ||
| 460 : | if (cl[1] != 't' || *diag_P(x) != 'N') { | ||
| 461 : | /* "trivially fast" when not triangular (<==> no 'diag' slot), | ||
| 462 : | or already *unit* triangular */ | ||
| 463 : | return (x); | ||
| 464 : | } | ||
| 465 : | else { /* triangular with diag='N'): now drop the diagonal */ | ||
| 466 : | /* duplicate, since chx will be modified: */ | ||
| 467 : | mmaechler | 2223 | CHM_SP chx = AS_CHM_SP__(duplicate(x)); |
| 468 : | maechler | 2125 | int uploT = (*uplo_P(x) == 'U') ? 1 : -1, |
| 469 : | Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; | ||
| 470 : | R_CheckStack(); | ||
| 471 : | |||
| 472 : | chm_diagN2U(chx, uploT, /* do_realloc */ FALSE); | ||
| 473 : | |||
| 474 : | return chm_sparse_to_SEXP(chx, /*dofree*/ 0/* or 1 ?? */, | ||
| 475 : | uploT, Rkind, "U", | ||
| 476 : | GET_SLOT(x, Matrix_DimNamesSym)); | ||
| 477 : | } | ||
| 478 : | } | ||
| 479 : | |||
| 480 : | bates | 1366 | SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j) |
| 481 : | { | ||
| 482 : | mmaechler | 2223 | CHM_SP chx = AS_CHM_SP__(x); |
| 483 : | bates | 1366 | int rsize = (isNull(i)) ? -1 : LENGTH(i), |
| 484 : | csize = (isNull(j)) ? -1 : LENGTH(j); | ||
| 485 : | maechler | 1736 | int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; |
| 486 : | maechler | 1960 | R_CheckStack(); |
| 487 : | bates | 1366 | |
| 488 : | if (rsize >= 0 && !isInteger(i)) | ||
| 489 : | error(_("Index i must be NULL or integer")); | ||
| 490 : | if (csize >= 0 && !isInteger(j)) | ||
| 491 : | error(_("Index j must be NULL or integer")); | ||
| 492 : | maechler | 1736 | |
| 493 : | dmbates | 2299 | return chm_sparse_to_SEXP(cholmod_l_submatrix(chx, INTEGER(i), rsize, |
| 494 : | maechler | 1375 | INTEGER(j), csize, |
| 495 : | bates | 1366 | TRUE, TRUE, &c), |
| 496 : | maechler | 1736 | 1, 0, Rkind, "", |
| 497 : | /* FIXME: drops dimnames */ R_NilValue); | ||
| 498 : | bates | 1366 | } |
| 499 : | bates | 2049 | |
| 500 : | SEXP Csparse_MatrixMarket(SEXP x, SEXP fname) | ||
| 501 : | { | ||
| 502 : | FILE *f = fopen(CHAR(asChar(fname)), "w"); | ||
| 503 : | |||
| 504 : | if (!f) | ||
| 505 : | error(_("failure to open file \"%s\" for writing"), | ||
| 506 : | CHAR(asChar(fname))); | ||
| 507 : | dmbates | 2299 | if (!cholmod_l_write_sparse(f, AS_CHM_SP(x), |
| 508 : | maechler | 2120 | (CHM_SP)NULL, (char*) NULL, &c)) |
| 509 : | dmbates | 2299 | error(_("cholmod_l_write_sparse returned error code")); |
| 510 : | bates | 2049 | fclose(f); |
| 511 : | return R_NilValue; | ||
| 512 : | } | ||
| 513 : | maechler | 2137 | |
| 514 : | |||
| 515 : | /** | ||
| 516 : | * Extract the diagonal entries from *triangular* Csparse matrix __or__ a | ||
| 517 : | * cholmod_sparse factor (LDL = TRUE). | ||
| 518 : | * | ||
| 519 : | * @param n dimension of the matrix. | ||
| 520 : | * @param x_p 'p' (column pointer) slot contents | ||
| 521 : | * @param x_x 'x' (non-zero entries) slot contents | ||
| 522 : | mmaechler | 2175 | * @param perm 'perm' (= permutation vector) slot contents; only used for "diagBack" |
| 523 : | maechler | 2137 | * @param resultKind a (SEXP) string indicating which kind of result is desired. |
| 524 : | * | ||
| 525 : | * @return a SEXP, either a (double) number or a length n-vector of diagonal entries | ||
| 526 : | */ | ||
| 527 : | SEXP diag_tC_ptr(int n, int *x_p, double *x_x, int *perm, SEXP resultKind) | ||
| 528 : | /* ^^^^^^ FIXME[Generalize] to int / ... */ | ||
| 529 : | { | ||
| 530 : | const char* res_ch = CHAR(STRING_ELT(resultKind,0)); | ||
| 531 : | enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log | ||
| 532 : | } res_kind = ((!strcmp(res_ch, "trace")) ? trace : | ||
| 533 : | ((!strcmp(res_ch, "sumLog")) ? sum_log : | ||
| 534 : | ((!strcmp(res_ch, "prod")) ? prod : | ||
| 535 : | ((!strcmp(res_ch, "diag")) ? diag : | ||
| 536 : | ((!strcmp(res_ch, "diagBack")) ? diag_backpermuted : | ||
| 537 : | -1))))); | ||
| 538 : | int i, n_x, i_from = 0; | ||
| 539 : | SEXP ans = PROTECT(allocVector(REALSXP, | ||
| 540 : | /* ^^^^ FIXME[Generalize] */ | ||
| 541 : | (res_kind == diag || | ||
| 542 : | res_kind == diag_backpermuted) ? n : 1)); | ||
| 543 : | double *v = REAL(ans); | ||
| 544 : | /* ^^^^^^ ^^^^ FIXME[Generalize] */ | ||
| 545 : | |||
| 546 : | #define for_DIAG(v_ASSIGN) \ | ||
| 547 : | for(i = 0; i < n; i++, i_from += n_x) { \ | ||
| 548 : | /* looking at i-th column */ \ | ||
| 549 : | n_x = x_p[i+1] - x_p[i];/* #{entries} in this column */ \ | ||
| 550 : | v_ASSIGN; \ | ||
| 551 : | } | ||
| 552 : | |||
| 553 : | /* NOTA BENE: we assume -- uplo = "L" i.e. lower triangular matrix | ||
| 554 : | * for uplo = "U" (makes sense with a "dtCMatrix" !), | ||
| 555 : | * should use x_x[i_from + (nx - 1)] instead of x_x[i_from], | ||
| 556 : | * where nx = (x_p[i+1] - x_p[i]) | ||
| 557 : | */ | ||
| 558 : | |||
| 559 : | switch(res_kind) { | ||
| 560 : | case trace: | ||
| 561 : | v[0] = 0.; | ||
| 562 : | for_DIAG(v[0] += x_x[i_from]); | ||
| 563 : | break; | ||
| 564 : | |||
| 565 : | case sum_log: | ||
| 566 : | v[0] = 0.; | ||
| 567 : | for_DIAG(v[0] += log(x_x[i_from])); | ||
| 568 : | break; | ||
| 569 : | |||
| 570 : | case prod: | ||
| 571 : | v[0] = 1.; | ||
| 572 : | for_DIAG(v[0] *= x_x[i_from]); | ||
| 573 : | break; | ||
| 574 : | |||
| 575 : | case diag: | ||
| 576 : | for_DIAG(v[i] = x_x[i_from]); | ||
| 577 : | break; | ||
| 578 : | |||
| 579 : | case diag_backpermuted: | ||
| 580 : | for_DIAG(v[i] = x_x[i_from]); | ||
| 581 : | |||
| 582 : | mmaechler | 2175 | warning(_("resultKind = 'diagBack' (back-permuted) is experimental")); |
| 583 : | maechler | 2144 | /* now back_permute : */ |
| 584 : | for(i = 0; i < n; i++) { | ||
| 585 : | double tmp = v[i]; v[i] = v[perm[i]]; v[perm[i]] = tmp; | ||
| 586 : | /*^^^^ FIXME[Generalize] */ | ||
| 587 : | } | ||
| 588 : | maechler | 2137 | break; |
| 589 : | |||
| 590 : | default: /* -1 from above */ | ||
| 591 : | error("diag_tC(): invalid 'resultKind'"); | ||
| 592 : | /* Wall: */ ans = R_NilValue; v = REAL(ans); | ||
| 593 : | } | ||
| 594 : | |||
| 595 : | UNPROTECT(1); | ||
| 596 : | return ans; | ||
| 597 : | } | ||
| 598 : | |||
| 599 : | /** | ||
| 600 : | * Extract the diagonal entries from *triangular* Csparse matrix __or__ a | ||
| 601 : | * cholmod_sparse factor (LDL = TRUE). | ||
| 602 : | * | ||
| 603 : | * @param pslot 'p' (column pointer) slot of Csparse matrix/factor | ||
| 604 : | * @param xslot 'x' (non-zero entries) slot of Csparse matrix/factor | ||
| 605 : | mmaechler | 2175 | * @param perm_slot 'perm' (= permutation vector) slot of corresponding CHMfactor; |
| 606 : | * only used for "diagBack" | ||
| 607 : | maechler | 2137 | * @param resultKind a (SEXP) string indicating which kind of result is desired. |
| 608 : | * | ||
| 609 : | * @return a SEXP, either a (double) number or a length n-vector of diagonal entries | ||
| 610 : | */ | ||
| 611 : | SEXP diag_tC(SEXP pslot, SEXP xslot, SEXP perm_slot, SEXP resultKind) | ||
| 612 : | { | ||
| 613 : | int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */ | ||
| 614 : | *x_p = INTEGER(pslot), | ||
| 615 : | *perm = INTEGER(perm_slot); | ||
| 616 : | double *x_x = REAL(xslot); | ||
| 617 : | /* ^^^^^^ ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/ | ||
| 618 : | |||
| 619 : | return diag_tC_ptr(n, x_p, x_x, perm, resultKind); | ||
| 620 : | } |
| R-Forge@R-project.org | ViewVC Help |
| Powered by ViewVC 1.0.0 |
