The R Project SVN R-packages

#include "dsCMatrix.h"

SEXP dsCMatrix_validate(SEXP obj)
{
    SEXP val = symmetricMatrix_validate(obj);
    if(isString(val))
    return(val);
    else {
    /* FIXME needed? dsC* inherits from dgC* which does this in validate*/
    csc_check_column_sorting(obj);
    return ScalarLogical(1);
    }
}

SEXP dsCMatrix_chol(SEXP x, SEXP pivot)
{
    cholmod_factor
    *N = as_cholmod_factor(dsCMatrix_Cholesky(x, pivot,
                          ScalarLogical(FALSE),
                          ScalarLogical(FALSE)));
    /* Must use a copy; cholmod_factor_to_sparse modifies first arg. */
    cholmod_factor *Ncp = cholmod_copy_factor(N, &c);
    cholmod_sparse *L, *R;
    SEXP ans;
    
    L = cholmod_factor_to_sparse(Ncp, &c); cholmod_free_factor(&Ncp, &c);
    R = cholmod_transpose(L, /*values*/ 1, &c); cholmod_free_sparse(&L, &c);
    ans = PROTECT(chm_sparse_to_SEXP(R, /*cholmod_free*/ 1,
                     /*uploT*/ 1, /*diag*/ "N",
                     GET_SLOT(x, Matrix_DimNamesSym)));
    if (asLogical(pivot)) {
    SEXP piv = PROTECT(allocVector(INTSXP, N->n));
    int *dest = INTEGER(piv), *src = (int*)N->Perm, i;
    
    for (i = 0; i < N->n; i++) dest[i] = src[i] + 1;
    setAttrib(ans, install("pivot"), piv);
    /* FIXME: Because of the cholmod_factor -> S4 obj ->
     * cholmod_factor conversions, the value of N->minor will
     * always be N->n.  Change as_cholmod_factor and
     * chm_factor_as_SEXP to keep track of Minor.
     */
    setAttrib(ans, install("rank"), ScalarInteger((size_t) N->minor));
    UNPROTECT(1);
    }
    Free(N);
    UNPROTECT(1);
    return ans;
}

SEXP dsCMatrix_Cholesky(SEXP Ap, SEXP permP, SEXP LDLp, SEXP superP)
{
    char *fname = strdup("spdCholesky"); /* template for factorization name */
    int perm = asLogical(permP),
    LDL = asLogical(LDLp),
    super = asLogical(superP);
    SEXP Chol;
    cholmod_sparse *A;
    cholmod_factor *L;
    int sup, ll;

    if (super) fname[0] = 'S';
    if (perm) fname[1] = 'P';
    if (LDL) fname[2] = 'D';
    Chol = get_factors(Ap, "fname");
    if (Chol != R_NilValue) return Chol;
    A = as_cholmod_sparse(Ap);
    sup = c.supernodal;
    ll = c.final_ll;
    
    if (!A->stype) error("Non-symmetric matrix passed to dsCMatrix_chol");
    
    c.final_ll = !LDL;  /* leave as LL' or form LDL' */
    c.supernodal = super ? CHOLMOD_SUPERNODAL : CHOLMOD_SIMPLICIAL; 

    if (perm) {
    L = cholmod_analyze(A, &c); /* get fill-reducing permutation */
    } else {                /* require identity permutation */
    int nmethods = c.nmethods, ord0 = c.method[0].ordering,
        postorder = c.postorder;
    c.nmethods = 1;
    c.method[0].ordering = CHOLMOD_NATURAL;
    c.postorder = FALSE;
    L = cholmod_analyze(A, &c);
    c.nmethods = nmethods; c.method[0].ordering = ord0;
    c.postorder = postorder;
    }
    c.supernodal = sup; /* restore previous setting */
    c.final_ll = ll;
    if (!cholmod_factorize(A, L, &c))
    error(_("Cholesky factorization failed"));
    Free(A);
    Chol = set_factors(Ap, chm_factor_to_SEXP(L, 1), fname);
    free(fname);        /* note, this must be free, not Free */
    return Chol;
}

static
SEXP get_factor_pattern(SEXP obj, char *pat, int offset)
{
    SEXP facs = GET_SLOT(obj, Matrix_factorSym), nms;
    int i;

    /* Why should this be nessary?  Shouldn't nms have length 0 if facs does? */
    if (!LENGTH(facs)) return R_NilValue;
    nms = getAttrib(facs, R_NamesSymbol);
    for (i = 0; i < LENGTH(nms); i++) {
    char *nm = CHAR(STRING_ELT(nms, i));
    if (strlen(nm) > offset && !strcmp(pat + offset, nm + offset))
        return VECTOR_ELT(facs, i);
    }
    return R_NilValue;
}

SEXP dsCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
{
    SEXP Chol = get_factor_pattern(a, "spdCholesky", 3);
    cholmod_factor *L;
    cholmod_dense *cb = as_cholmod_dense(b), *cx;
    
    if (Chol == R_NilValue)
    Chol = dsCMatrix_Cholesky(a,
                  ScalarLogical(1),  /* permuted */
                  ScalarLogical(1),  /* LDL' */
                  ScalarLogical(0)); /* simplicial */
    L = as_cholmod_factor(Chol);
    cx = cholmod_solve(CHOLMOD_A, L, cb, &c);
    Free(cb); Free(L);
    return chm_dense_to_SEXP(cx, 1);
}

/* TODO: still needed with Csparse_to_Tsparse()  [ which does dsC -> dsT ] */
SEXP dsCMatrix_to_dgTMatrix(SEXP x)
{
    SEXP
    ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix"))),
    islot = GET_SLOT(x, Matrix_iSym),
    pslot = GET_SLOT(x, Matrix_pSym);
    int *ai, *aj, *iv = INTEGER(islot),
    j, jj, nnz = length(islot), nout,
    n = length(pslot) - 1,
    *p = INTEGER(pslot), pos;
    double *ax, *xv = REAL(GET_SLOT(x, Matrix_xSym));

    /* increment output count by number of off-diagonals */
    nout = nnz;
    for (j = 0; j < n; j++) {
    int p2 = p[j+1];
    for (jj = p[j]; jj < p2; jj++) {
        if (iv[jj] != j) nout++;
    }
    }
    SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
    SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nout));
    ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
    SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, nout));
    aj = INTEGER(GET_SLOT(ans, Matrix_jSym));
    SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nout));
    ax = REAL(GET_SLOT(ans, Matrix_xSym));
    pos = 0;
    for (j = 0; j < n; j++) {
    int p2 = p[j+1];
    for (jj = p[j]; jj < p2; jj++) {
        int ii = iv[jj];
        double xx = xv[jj];

        ai[pos] = ii; aj[pos] = j; ax[pos] = xx; pos++;
        if (ii != j) {
        aj[pos] = ii; ai[pos] = j; ax[pos] = xx; pos++;
        }
    }
    }
    UNPROTECT(1);
    return ans;
}