--- pkg/Matrix/src/Csparse.c 2011/08/05 19:52:10 2685 +++ pkg/Matrix/src/Csparse.c 2014/10/06 17:00:48 3011 @@ -37,67 +37,23 @@ return Csparse_validate_(x, FALSE); } + +#define _t_Csparse_validate +#include "t_Csparse_validate.c" + +#define _t_Csparse_sort +#include "t_Csparse_validate.c" + +// R: .validateCsparse(x, sort.if.needed = FALSE) : SEXP Csparse_validate2(SEXP x, SEXP maybe_modify) { return Csparse_validate_(x, asLogical(maybe_modify)); } -SEXP Csparse_validate_(SEXP x, Rboolean maybe_modify) -{ - /* NB: we do *NOT* check a potential 'x' slot here, at all */ - SEXP pslot = GET_SLOT(x, Matrix_pSym), - islot = GET_SLOT(x, Matrix_iSym); - Rboolean sorted, strictly; - int j, k, - *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), - nrow = dims[0], - ncol = dims[1], - *xp = INTEGER(pslot), - *xi = INTEGER(islot); - - if (length(pslot) != dims[1] + 1) - return mkString(_("slot p must have length = ncol(.) + 1")); - if (xp[0] != 0) - return mkString(_("first element of slot p must be zero")); - if (length(islot) < xp[ncol]) /* allow larger slots from over-allocation!*/ - return - mkString(_("last element of slot p must match length of slots i and x")); - for (j = 0; j < xp[ncol]; j++) { - if (xi[j] < 0 || xi[j] >= nrow) - return mkString(_("all row indices must be between 0 and nrow-1")); - } - sorted = TRUE; strictly = TRUE; - for (j = 0; j < ncol; j++) { - if (xp[j] > xp[j + 1]) - return mkString(_("slot p must be non-decreasing")); - if(sorted) /* only act if >= 2 entries in column j : */ - for (k = xp[j] + 1; k < xp[j + 1]; k++) { - if (xi[k] < xi[k - 1]) - sorted = FALSE; - else if (xi[k] == xi[k - 1]) - strictly = FALSE; - } - } - if (!sorted) { - if(maybe_modify) { - CHM_SP chx = (CHM_SP) alloca(sizeof(cholmod_sparse)); - R_CheckStack(); - as_cholmod_sparse(chx, x, FALSE, TRUE);/*-> cholmod_l_sort() ! */ - /* as chx = AS_CHM_SP__(x) but ^^^^ sorting x in_place !!! */ - - /* Now re-check that row indices are *strictly* increasing - * (and not just increasing) within each column : */ - for (j = 0; j < ncol; j++) { - for (k = xp[j] + 1; k < xp[j + 1]; k++) - if (xi[k] == xi[k - 1]) - return mkString(_("slot i is not *strictly* increasing inside a column (even after cholmod_l_sort)")); - } - } else { /* no modifying sorting : */ - return mkString(_("row indices are not sorted within columns")); - } - } else if(!strictly) { /* sorted, but not strictly */ - return mkString(_("slot i is not *strictly* increasing inside a column")); - } - return ScalarLogical(1); +// R: Matrix:::.sortCsparse(x) : +SEXP Csparse_sort (SEXP x) { + int ok = Csparse_sort_2(x, TRUE); // modifying x directly + if(!ok) warning(_("Csparse_sort(x): x is not a valid (apart from sorting) CsparseMatrix")); + return x; } SEXP Rsparse_validate(SEXP x) @@ -189,7 +145,8 @@ if(cl_x[2] != 'C') error(_("not a CsparseMatrix")); int nnz = LENGTH(GET_SLOT(x, Matrix_iSym)); SEXP ans; - char *ncl = strdup(cl_x); + char *ncl = alloca(strlen(cl_x) + 1); /* not much memory required */ + strcpy(ncl, cl_x); double *dx_x; int *ix_x; ncl[0] = (r_kind == x_double ? 'd' : (r_kind == x_logical ? 'l' : @@ -229,11 +186,15 @@ return ans; } -SEXP Csparse_to_matrix(SEXP x) +SEXP Csparse_to_matrix(SEXP x, SEXP chk) { - return chm_dense_to_matrix(cholmod_sparse_to_dense(AS_CHM_SP__(x), &c), + return chm_dense_to_matrix(cholmod_sparse_to_dense(AS_CHM_SP2(x, asLogical(chk)), &c), 1 /*do_free*/, GET_SLOT(x, Matrix_DimNamesSym)); } +SEXP Csparse_to_vector(SEXP x) +{ + return chm_dense_to_vector(cholmod_sparse_to_dense(AS_CHM_SP__(x), &c), 1); +} SEXP Csparse_to_Tsparse(SEXP x, SEXP tri) { @@ -266,11 +227,15 @@ SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo) { + int *adims = INTEGER(GET_SLOT(x, Matrix_DimSym)), n = adims[0]; + if(n != adims[1]) { + error(_("Csparse_general_to_symmetric(): matrix is not square!")); + return R_NilValue; /* -Wall */ + } CHM_SP chx = AS_CHM_SP__(x), chgx; int uploT = (*CHAR(STRING_ELT(uplo,0)) == 'U') ? 1 : -1; int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; R_CheckStack(); - chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c); /* xtype: pattern, "real", complex or .. */ return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", @@ -488,7 +453,7 @@ if(!cholmod_drop(dtol, ans, &c)) error(_("cholmod_drop() failed")); - return chm_sparse_to_SEXP(ans, 1, + return chm_sparse_to_SEXP(ans, 1, tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0, Rkind, tr ? diag_P(x) : "", GET_SLOT(x, Matrix_DimNamesSym)); @@ -566,16 +531,19 @@ } else { /* triangular with diag='N'): now drop the diagonal */ /* duplicate, since chx will be modified: */ - CHM_SP chx = AS_CHM_SP__(duplicate(x)); + SEXP xx = PROTECT(duplicate(x)); + CHM_SP chx = AS_CHM_SP__(xx); int uploT = (*uplo_P(x) == 'U') ? 1 : -1, Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0; R_CheckStack(); chm_diagN2U(chx, uploT, /* do_realloc */ FALSE); - return chm_sparse_to_SEXP(chx, /*dofree*/ 0/* or 1 ?? */, - uploT, Rkind, "U", - GET_SLOT(x, Matrix_DimNamesSym)); + SEXP ans = chm_sparse_to_SEXP(chx, /*dofree*/ 0/* or 1 ?? */, + uploT, Rkind, "U", + GET_SLOT(x, Matrix_DimNamesSym)); + UNPROTECT(1);// only now ! + return ans; } } @@ -601,16 +569,23 @@ if (csize >= 0 && !isInteger(j)) error(_("Index j must be NULL or integer")); - if (chx->stype) /* symmetricMatrix */ +#define CHM_SUB(_M_, _i_, _j_) \ + cholmod_submatrix(_M_, \ + (rsize < 0) ? NULL : INTEGER(_i_), rsize, \ + (csize < 0) ? NULL : INTEGER(_j_), csize, \ + TRUE, TRUE, &c) + CHM_SP ans; + if (!chx->stype) {/* non-symmetric Matrix */ + ans = CHM_SUB(chx, i, j); + } + else { /* for now, cholmod_submatrix() only accepts "generalMatrix" */ - chx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c); - - return chm_sparse_to_SEXP(cholmod_submatrix(chx, - (rsize < 0) ? NULL : INTEGER(i), rsize, - (csize < 0) ? NULL : INTEGER(j), csize, - TRUE, TRUE, &c), - 1, 0, Rkind, "", - /* FIXME: drops dimnames */ R_NilValue); + CHM_SP tmp = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c); + ans = CHM_SUB(tmp, i, j); + cholmod_free_sparse(&tmp, &c); + } + // "FIXME": currently dropping dimnames, and adding them afterwards in R : + return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", /* dimnames: */ R_NilValue); } #define _d_Csp_ @@ -657,40 +632,49 @@ * * @return a SEXP, either a (double) number or a length n-vector of diagonal entries */ -SEXP diag_tC_ptr(int n, int *x_p, double *x_x, int *perm, SEXP resultKind) +SEXP diag_tC_ptr(int n, int *x_p, double *x_x, Rboolean is_U, int *perm, /* ^^^^^^ FIXME[Generalize] to int / ... */ + SEXP resultKind) { const char* res_ch = CHAR(STRING_ELT(resultKind,0)); - enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log + enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log, min, max, range } res_kind = ((!strcmp(res_ch, "trace")) ? trace : ((!strcmp(res_ch, "sumLog")) ? sum_log : ((!strcmp(res_ch, "prod")) ? prod : - ((!strcmp(res_ch, "diag")) ? diag : - ((!strcmp(res_ch, "diagBack")) ? diag_backpermuted : - -1))))); - int i, n_x, i_from = 0; + ((!strcmp(res_ch, "min")) ? min : + ((!strcmp(res_ch, "max")) ? max : + ((!strcmp(res_ch, "range")) ? range : + ((!strcmp(res_ch, "diag")) ? diag : + ((!strcmp(res_ch, "diagBack")) ? diag_backpermuted : + -1)))))))); + int i, n_x, i_from; SEXP ans = PROTECT(allocVector(REALSXP, /* ^^^^ FIXME[Generalize] */ (res_kind == diag || - res_kind == diag_backpermuted) ? n : 1)); + res_kind == diag_backpermuted) ? n : + (res_kind == range ? 2 : 1))); double *v = REAL(ans); /* ^^^^^^ ^^^^ FIXME[Generalize] */ -#define for_DIAG(v_ASSIGN) \ - for(i = 0; i < n; i++, i_from += n_x) { \ - /* looking at i-th column */ \ + i_from = (is_U ? -1 : 0); + +#define for_DIAG(v_ASSIGN) \ + for(i = 0; i < n; i++) { \ + /* looking at i-th column */ \ n_x = x_p[i+1] - x_p[i];/* #{entries} in this column */ \ - v_ASSIGN; \ + if( is_U) i_from += n_x; \ + v_ASSIGN; \ + if(!is_U) i_from += n_x; \ } /* NOTA BENE: we assume -- uplo = "L" i.e. lower triangular matrix * for uplo = "U" (makes sense with a "dtCMatrix" !), - * should use x_x[i_from + (nx - 1)] instead of x_x[i_from], - * where nx = (x_p[i+1] - x_p[i]) + * should use x_x[i_from + (n_x - 1)] instead of x_x[i_from], + * where n_x = (x_p[i+1] - x_p[i]) */ switch(res_kind) { - case trace: + case trace: // = sum v[0] = 0.; for_DIAG(v[0] += x_x[i_from]); break; @@ -705,6 +689,23 @@ for_DIAG(v[0] *= x_x[i_from]); break; + case min: + v[0] = R_PosInf; + for_DIAG(if(v[0] > x_x[i_from]) v[0] = x_x[i_from]); + break; + + case max: + v[0] = R_NegInf; + for_DIAG(if(v[0] < x_x[i_from]) v[0] = x_x[i_from]); + break; + + case range: + v[0] = R_PosInf; + v[1] = R_NegInf; + for_DIAG(if(v[0] > x_x[i_from]) v[0] = x_x[i_from]; + if(v[1] < x_x[i_from]) v[1] = x_x[i_from]); + break; + case diag: for_DIAG(v[i] = x_x[i_from]); break; @@ -712,7 +713,8 @@ case diag_backpermuted: for_DIAG(v[i] = x_x[i_from]); - warning(_("resultKind = 'diagBack' (back-permuted) is experimental")); + warning(_("%s = '%s' (back-permuted) is experimental"), + "resultKind", "diagBack"); /* now back_permute : */ for(i = 0; i < n; i++) { double tmp = v[i]; v[i] = v[perm[i]]; v[perm[i]] = tmp; @@ -733,6 +735,7 @@ * Extract the diagonal entries from *triangular* Csparse matrix __or__ a * cholmod_sparse factor (LDL = TRUE). * + * @param obj -- now a cholmod_sparse factor or a dtCMatrix * @param pslot 'p' (column pointer) slot of Csparse matrix/factor * @param xslot 'x' (non-zero entries) slot of Csparse matrix/factor * @param perm_slot 'perm' (= permutation vector) slot of corresponding CHMfactor; @@ -741,17 +744,27 @@ * * @return a SEXP, either a (double) number or a length n-vector of diagonal entries */ -SEXP diag_tC(SEXP pslot, SEXP xslot, SEXP perm_slot, SEXP resultKind) +SEXP diag_tC(SEXP obj, SEXP resultKind) { + + SEXP + pslot = GET_SLOT(obj, Matrix_pSym), + xslot = GET_SLOT(obj, Matrix_xSym); + Rboolean is_U = (R_has_slot(obj, Matrix_uploSym) && + *CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))) == 'U'); int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */ - *x_p = INTEGER(pslot), - *perm = INTEGER(perm_slot); + *x_p = INTEGER(pslot), pp = -1, *perm; double *x_x = REAL(xslot); /* ^^^^^^ ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/ - return diag_tC_ptr(n, x_p, x_x, perm, resultKind); + if(R_has_slot(obj, Matrix_permSym)) + perm = INTEGER(GET_SLOT(obj, Matrix_permSym)); + else perm = &pp; + + return diag_tC_ptr(n, x_p, x_x, is_U, perm, resultKind); } + /** * Create a Csparse matrix object from indices and/or pointers. *