The R Project SVN R-packages

#include "Pattern.h"

cholmod_sparse *factor_to_pattern(SEXP fact)

{

    int *vals = INTEGER(fact), i, n = LENGTH(fact),

	nlev = LENGTH(getAttrib(fact, R_LevelsSymbol));

    cholmod_sparse *A;

    if (!isFactor(fact))

	error("`fact' must be a factor");

    A = cholmod_allocate_sparse((size_t) nlev, (size_t) n, (size_t) n,

				TRUE, TRUE, 0, CHOLMOD_PATTERN, &c);

    for (i = 0; i <= n; i++) ((int *) (A->p))[i] = i;

    for (i = 0; i < n; i++) ((int *) (A->i))[i] = vals[i] - 1;

    if (cholmod_check_sparse(A, &c) < 0)

	error("sparse Pattern matrix has invalid structure");

    return A;

}

SEXP factor_prod(SEXP f1, SEXP f2)

{

    SEXP ans = PROTECT(allocVector(VECSXP, 3));

    int *dims, n = LENGTH(f1), i, nl1, nl2, nmax, super;

    cholmod_triplet *A;

    cholmod_sparse *B;

    cholmod_factor *F;

    if (!isFactor(f1) || !isFactor(f2) || LENGTH(f2) != n)

	error("f1 and f2 must be factors of the same length");

    nl1 = LENGTH(getAttrib(f1, R_LevelsSymbol));

    nl2 = LENGTH(getAttrib(f2, R_LevelsSymbol));

    A = cholmod_allocate_triplet((size_t) nl1, (size_t) nl2, (size_t) n,

				 0, CHOLMOD_PATTERN, &c);

    for (i = 0; i < n; i++) {

	((int *)(A->i))[i] = INTEGER(f1)[i] - 1;

	((int *)(A->j))[i] = INTEGER(f2)[i] - 1;

    }

    A->nnz = n;

    nmax = nl1 * nl2;

    if (n < nmax) nmax = n;

    B = cholmod_triplet_to_sparse(A, nmax, &c);

    cholmod_free_triplet(&A, &c);

				/*  force a simplicial factorization */

    super = c.supernodal;

    c.supernodal = CHOLMOD_SIMPLICIAL;

    F = cholmod_analyze(B, &c); 

    c.supernodal = super;

    SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, 2));

    dims = INTEGER(VECTOR_ELT(ans, 0));

    dims[0] = nl1; dims[1] = nl2;

    SET_VECTOR_ELT(ans, 1, allocVector(INTSXP, nl1));

    Memcpy(INTEGER(VECTOR_ELT(ans, 1)), (int *) F->Perm, nl1);

    SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, nl1));

    Memcpy(INTEGER(VECTOR_ELT(ans, 2)), (int *) F->ColCount, nl1);

    Rprintf("Ordering used: %d\n", F->ordering);

    Rprintf("is_super: %d\n", F->is_super);

    cholmod_free_sparse(&B, &c);

    cholmod_free_factor(&F, &c);

    UNPROTECT(1);

    return ans;

}