SCM Repository
Annotation of /pkg/src/Csparse.c
Parent Directory
| 
Revision Log
Revision 1461 - (view) (download) (as text)
| 1 : | bates | 1218 | /* Sparse matrices in compressed column-oriented form */ |
| 2 : | bates | 922 | #include "Csparse.h" |
| 3 : | #include "chm_common.h" | ||
| 4 : | |||
| 5 : | SEXP Csparse_validate(SEXP x) | ||
| 6 : | { | ||
| 7 : | SEXP pslot = GET_SLOT(x, Matrix_pSym), | ||
| 8 : | islot = GET_SLOT(x, Matrix_iSym); | ||
| 9 : | int j, ncol = length(pslot) - 1, | ||
| 10 : | *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), | ||
| 11 : | nrow, *xp = INTEGER(pslot), | ||
| 12 : | *xi = INTEGER(islot); | ||
| 13 : | |||
| 14 : | nrow = dims[0]; | ||
| 15 : | if (length(pslot) <= 0) | ||
| 16 : | return mkString(_("slot p must have length > 0")); | ||
| 17 : | if (xp[0] != 0) | ||
| 18 : | return mkString(_("first element of slot p must be zero")); | ||
| 19 : | if (length(islot) != xp[ncol]) | ||
| 20 : | return mkString(_("last element of slot p must match length of slots i and x")); | ||
| 21 : | for (j = 0; j < ncol; j++) { | ||
| 22 : | if (xp[j] > xp[j+1]) | ||
| 23 : | return mkString(_("slot p must be non-decreasing")); | ||
| 24 : | } | ||
| 25 : | for (j = 0; j < length(islot); j++) { | ||
| 26 : | if (xi[j] < 0 || xi[j] >= nrow) | ||
| 27 : | return mkString(_("all row indices must be between 0 and nrow-1")); | ||
| 28 : | } | ||
| 29 : | return ScalarLogical(1); | ||
| 30 : | } | ||
| 31 : | |||
| 32 : | bates | 1059 | SEXP Csparse_to_dense(SEXP x) |
| 33 : | { | ||
| 34 : | cholmod_sparse *chxs = as_cholmod_sparse(x); | ||
| 35 : | cholmod_dense *chxd = cholmod_sparse_to_dense(chxs, &c); | ||
| 36 : | |||
| 37 : | bates | 1141 | Free(chxs); |
| 38 : | bates | 1059 | return chm_dense_to_SEXP(chxd, 1); |
| 39 : | } | ||
| 40 : | |||
| 41 : | bates | 1371 | SEXP Csparse_to_logical(SEXP x, SEXP tri) |
| 42 : | { | ||
| 43 : | cholmod_sparse *chxs = as_cholmod_sparse(x); | ||
| 44 : | cholmod_sparse | ||
| 45 : | *chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c); | ||
| 46 : | int uploT = 0; char *diag = ""; | ||
| 47 : | |||
| 48 : | Free(chxs); | ||
| 49 : | if (asLogical(tri)) { /* triangular sparse matrices */ | ||
| 50 : | uploT = (strcmp(CHAR(asChar(GET_SLOT(x, Matrix_uploSym))), "U")) ? | ||
| 51 : | -1 : 1; | ||
| 52 : | diag = CHAR(asChar(GET_SLOT(x, Matrix_diagSym))); | ||
| 53 : | } | ||
| 54 : | return chm_sparse_to_SEXP(chxcp, 1, uploT, diag, | ||
| 55 : | GET_SLOT(x, Matrix_DimNamesSym)); | ||
| 56 : | } | ||
| 57 : | |||
| 58 : | bates | 1366 | SEXP Csparse_to_matrix(SEXP x) |
| 59 : | bates | 922 | { |
| 60 : | maechler | 925 | cholmod_sparse *chxs = as_cholmod_sparse(x); |
| 61 : | bates | 1366 | cholmod_dense *chxd = cholmod_sparse_to_dense(chxs, &c); |
| 62 : | |||
| 63 : | Free(chxs); | ||
| 64 : | return chm_dense_to_matrix(chxd, 1, | ||
| 65 : | GET_SLOT(x, Matrix_DimNamesSym)); | ||
| 66 : | } | ||
| 67 : | |||
| 68 : | SEXP Csparse_to_Tsparse(SEXP x, SEXP tri) | ||
| 69 : | { | ||
| 70 : | cholmod_sparse *chxs = as_cholmod_sparse(x); | ||
| 71 : | bates | 922 | cholmod_triplet *chxt = cholmod_sparse_to_triplet(chxs, &c); |
| 72 : | bates | 1366 | int uploT = 0; char *diag = ""; |
| 73 : | bates | 922 | |
| 74 : | bates | 1141 | Free(chxs); |
| 75 : | bates | 1366 | if (asLogical(tri)) { /* triangular sparse matrices */ |
| 76 : | bates | 1419 | uploT = (*uplo_P(x) == 'U') ? -1 : 1; |
| 77 : | diag = diag_P(x); | ||
| 78 : | bates | 1366 | } |
| 79 : | return chm_triplet_to_SEXP(chxt, 1, uploT, diag, | ||
| 80 : | bates | 1371 | GET_SLOT(x, Matrix_DimNamesSym)); |
| 81 : | bates | 922 | } |
| 82 : | |||
| 83 : | bates | 1448 | /* this used to be called sCMatrix_to_gCMatrix(..) [in ./dsCMatrix.c ]: */ |
| 84 : | bates | 1371 | SEXP Csparse_symmetric_to_general(SEXP x) |
| 85 : | { | ||
| 86 : | cholmod_sparse *chx = as_cholmod_sparse(x), *chgx; | ||
| 87 : | |||
| 88 : | if (!(chx->stype)) | ||
| 89 : | error(_("Nonsymmetric matrix in Csparse_symmeteric_to_general")); | ||
| 90 : | maechler | 1375 | chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c); |
| 91 : | /* xtype: pattern, "real", complex or .. */ | ||
| 92 : | bates | 1371 | Free(chx); |
| 93 : | return chm_sparse_to_SEXP(chgx, 1, 0, "", | ||
| 94 : | GET_SLOT(x, Matrix_DimNamesSym)); | ||
| 95 : | } | ||
| 96 : | |||
| 97 : | bates | 1369 | SEXP Csparse_transpose(SEXP x, SEXP tri) |
| 98 : | bates | 922 | { |
| 99 : | cholmod_sparse *chx = as_cholmod_sparse(x); | ||
| 100 : | cholmod_sparse *chxt = cholmod_transpose(chx, (int) chx->xtype, &c); | ||
| 101 : | bates | 1366 | SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp; |
| 102 : | bates | 1369 | int uploT = 0; char *diag = ""; |
| 103 : | |||
| 104 : | bates | 1141 | Free(chx); |
| 105 : | bates | 1366 | tmp = VECTOR_ELT(dn, 0); /* swap the dimnames */ |
| 106 : | SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1)); | ||
| 107 : | SET_VECTOR_ELT(dn, 1, tmp); | ||
| 108 : | UNPROTECT(1); | ||
| 109 : | bates | 1369 | if (asLogical(tri)) { /* triangular sparse matrices */ |
| 110 : | bates | 1419 | uploT = (*uplo_P(x) == 'U') ? -1 : 1; |
| 111 : | diag = diag_P(x); | ||
| 112 : | bates | 1369 | } |
| 113 : | return chm_sparse_to_SEXP(chxt, 1, uploT, diag, dn); | ||
| 114 : | bates | 922 | } |
| 115 : | |||
| 116 : | SEXP Csparse_Csparse_prod(SEXP a, SEXP b) | ||
| 117 : | { | ||
| 118 : | bates | 1059 | cholmod_sparse *cha = as_cholmod_sparse(a), |
| 119 : | *chb = as_cholmod_sparse(b); | ||
| 120 : | cholmod_sparse *chc = cholmod_ssmult(cha, chb, 0, cha->xtype, 1, &c); | ||
| 121 : | bates | 1366 | SEXP dn = allocVector(VECSXP, 2); |
| 122 : | bates | 922 | |
| 123 : | bates | 1141 | Free(cha); Free(chb); |
| 124 : | bates | 1366 | SET_VECTOR_ELT(dn, 0, /* establish dimnames */ |
| 125 : | duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0))); | ||
| 126 : | SET_VECTOR_ELT(dn, 1, | ||
| 127 : | duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1))); | ||
| 128 : | return chm_sparse_to_SEXP(chc, 1, 0, "", dn); | ||
| 129 : | bates | 922 | } |
| 130 : | |||
| 131 : | SEXP Csparse_dense_prod(SEXP a, SEXP b) | ||
| 132 : | { | ||
| 133 : | cholmod_sparse *cha = as_cholmod_sparse(a); | ||
| 134 : | bates | 1461 | cholmod_dense *chb = as_cholmod_dense(PROTECT(mMatrix_as_dgeMatrix(b))); |
| 135 : | bates | 1448 | cholmod_dense *chc = |
| 136 : | cholmod_allocate_dense(cha->nrow, chb->ncol, cha->nrow, chb->xtype, &c); | ||
| 137 : | double alpha[] = {1,0}, beta[] = {0,0}; | ||
| 138 : | bates | 922 | |
| 139 : | bates | 1448 | cholmod_sdmult(cha, 0, alpha, beta, chb, chc, &c); |
| 140 : | bates | 1141 | Free(cha); Free(chb); |
| 141 : | bates | 1461 | UNPROTECT(1); |
| 142 : | bates | 923 | return chm_dense_to_SEXP(chc, 1); |
| 143 : | bates | 922 | } |
| 144 : | maechler | 925 | |
| 145 : | bates | 1067 | SEXP Csparse_dense_crossprod(SEXP a, SEXP b) |
| 146 : | { | ||
| 147 : | cholmod_sparse *cha = as_cholmod_sparse(a); | ||
| 148 : | bates | 1461 | cholmod_dense *chb = as_cholmod_dense(PROTECT(mMatrix_as_dgeMatrix(b))); |
| 149 : | bates | 1448 | cholmod_dense *chc = |
| 150 : | cholmod_allocate_dense(cha->ncol, chb->ncol, cha->ncol, chb->xtype, &c); | ||
| 151 : | double alpha[] = {1,0}, beta[] = {0,0}; | ||
| 152 : | bates | 1067 | |
| 153 : | bates | 1448 | cholmod_sdmult(cha, 1, alpha, beta, chb, chc, &c); |
| 154 : | bates | 1067 | Free(cha); Free(chb); |
| 155 : | bates | 1461 | UNPROTECT(1); |
| 156 : | bates | 1067 | return chm_dense_to_SEXP(chc, 1); |
| 157 : | } | ||
| 158 : | |||
| 159 : | bates | 928 | SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet) |
| 160 : | bates | 922 | { |
| 161 : | maechler | 957 | int trip = asLogical(triplet), |
| 162 : | tr = asLogical(trans); /* gets reversed because _aat is tcrossprod */ | ||
| 163 : | bates | 928 | cholmod_triplet |
| 164 : | *cht = trip ? as_cholmod_triplet(x) : (cholmod_triplet*) NULL; | ||
| 165 : | cholmod_sparse *chcp, *chxt, | ||
| 166 : | *chx = trip ? cholmod_triplet_to_sparse(cht, cht->nnz, &c) | ||
| 167 : | : as_cholmod_sparse(x); | ||
| 168 : | bates | 1366 | SEXP dn = PROTECT(allocVector(VECSXP, 2)); |
| 169 : | bates | 922 | |
| 170 : | bates | 923 | if (!tr) |
| 171 : | bates | 1353 | chxt = cholmod_transpose(chx, chx->xtype, &c); |
| 172 : | bates | 928 | chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c); |
| 173 : | maechler | 957 | if(!chcp) |
| 174 : | error("Csparse_crossprod(): error return from cholmod_aat()"); | ||
| 175 : | bates | 1360 | cholmod_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c); |
| 176 : | chcp->stype = 1; | ||
| 177 : | bates | 930 | if (trip) { |
| 178 : | cholmod_free_sparse(&chx, &c); | ||
| 179 : | bates | 1141 | Free(cht); |
| 180 : | bates | 930 | } else { |
| 181 : | bates | 1141 | Free(chx); |
| 182 : | bates | 930 | } |
| 183 : | bates | 923 | if (!tr) cholmod_free_sparse(&chxt, &c); |
| 184 : | bates | 1366 | /* create dimnames */ |
| 185 : | SET_VECTOR_ELT(dn, 0, | ||
| 186 : | duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), | ||
| 187 : | (tr) ? 1 : 0))); | ||
| 188 : | SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(dn, 0))); | ||
| 189 : | UNPROTECT(1); | ||
| 190 : | return chm_sparse_to_SEXP(chcp, 1, 0, "", dn); | ||
| 191 : | bates | 922 | } |
| 192 : | bates | 923 | |
| 193 : | bates | 1218 | SEXP Csparse_horzcat(SEXP x, SEXP y) |
| 194 : | { | ||
| 195 : | cholmod_sparse *chx = as_cholmod_sparse(x), | ||
| 196 : | *chy = as_cholmod_sparse(y), *ans; | ||
| 197 : | maechler | 1375 | |
| 198 : | bates | 1218 | ans = cholmod_horzcat(chx, chy, 1, &c); |
| 199 : | Free(chx); Free(chy); | ||
| 200 : | bates | 1366 | /* FIXME: currently drops dimnames */ |
| 201 : | return chm_sparse_to_SEXP(ans, 1, 0, "", R_NilValue); | ||
| 202 : | bates | 1218 | } |
| 203 : | |||
| 204 : | SEXP Csparse_vertcat(SEXP x, SEXP y) | ||
| 205 : | { | ||
| 206 : | cholmod_sparse *chx = as_cholmod_sparse(x), | ||
| 207 : | *chy = as_cholmod_sparse(y), *ans; | ||
| 208 : | maechler | 1375 | |
| 209 : | bates | 1218 | ans = cholmod_vertcat(chx, chy, 1, &c); |
| 210 : | Free(chx); Free(chy); | ||
| 211 : | bates | 1366 | /* FIXME: currently drops dimnames */ |
| 212 : | return chm_sparse_to_SEXP(ans, 1, 0, "", R_NilValue); | ||
| 213 : | bates | 1218 | } |
| 214 : | bates | 1265 | |
| 215 : | SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2) | ||
| 216 : | { | ||
| 217 : | cholmod_sparse *chx = as_cholmod_sparse(x), *ans; | ||
| 218 : | |||
| 219 : | ans = cholmod_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c); | ||
| 220 : | Free(chx); | ||
| 221 : | bates | 1366 | return chm_sparse_to_SEXP(ans, 1, 0, "", R_NilValue); |
| 222 : | bates | 1265 | } |
| 223 : | bates | 1366 | |
| 224 : | SEXP Csparse_diagU2N(SEXP x) | ||
| 225 : | { | ||
| 226 : | cholmod_sparse *chx = as_cholmod_sparse(x); | ||
| 227 : | cholmod_sparse *eye = cholmod_speye(chx->nrow, chx->ncol, chx->xtype, &c); | ||
| 228 : | double one[] = {1, 0}; | ||
| 229 : | cholmod_sparse *ans = cholmod_add(chx, eye, one, one, TRUE, TRUE, &c); | ||
| 230 : | int uploT = (strcmp(CHAR(asChar(GET_SLOT(x, Matrix_uploSym))), "U")) ? | ||
| 231 : | -1 : 1; | ||
| 232 : | |||
| 233 : | Free(chx); cholmod_free_sparse(&eye, &c); | ||
| 234 : | return chm_sparse_to_SEXP(ans, 1, uploT, "N", | ||
| 235 : | duplicate(GET_SLOT(x, Matrix_DimNamesSym))); | ||
| 236 : | } | ||
| 237 : | |||
| 238 : | SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j) | ||
| 239 : | { | ||
| 240 : | cholmod_sparse *chx = as_cholmod_sparse(x); | ||
| 241 : | int rsize = (isNull(i)) ? -1 : LENGTH(i), | ||
| 242 : | csize = (isNull(j)) ? -1 : LENGTH(j); | ||
| 243 : | |||
| 244 : | if (rsize >= 0 && !isInteger(i)) | ||
| 245 : | error(_("Index i must be NULL or integer")); | ||
| 246 : | if (csize >= 0 && !isInteger(j)) | ||
| 247 : | error(_("Index j must be NULL or integer")); | ||
| 248 : | return chm_sparse_to_SEXP(cholmod_submatrix(chx, INTEGER(i), rsize, | ||
| 249 : | maechler | 1375 | INTEGER(j), csize, |
| 250 : | bates | 1366 | TRUE, TRUE, &c), |
| 251 : | 1, 0, "", R_NilValue); | ||
| 252 : | } |
| root@r-forge.r-project.org | ViewVC Help |
| Powered by ViewVC 1.0.0 |
