The R Project SVN R-packages

                /* LDL' factorization of a logical,
                 * sparse column-oriented matrix */
#include "lCholCMatrix.h"

SEXP lCholCMatrix_validate(SEXP x)
{
    int i, n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];
    SEXP perm = GET_SLOT(x, Matrix_permSym),
    Parent = GET_SLOT(x, Matrix_ParentSym);

    if (length(perm) != n)
    return mkString(_("slot perm must have length n"));
    if (length(Parent) != n)
    return mkString(_("slot Parent must have length n"));
    if (!R_ldl_valid_perm(n, INTEGER(perm))) 
    return mkString(_("slot perm is not a valid 0-based permutation"));
    for (i = 0; i < n; i++) {
    int pari = INTEGER(Parent)[i];
    if (pari < -1 || pari > (n - 1))
        return mkString(_("an element of the Parent array is not in range [-1,n-1]"));
    }

    return ScalarLogical(1);
}

/** 
 * Create the structure of the inverse of L from the LDL' factorization.
 * 
 * @param x Pointer to a lCholCMatrix object
 * 
 * @return An ltCMatrix object representing L^{-1}
 */
SEXP lCholCMatrix_solve(SEXP x)
{
    SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("ltCMatrix")));
    SEXP Parent = GET_SLOT(x, Matrix_ParentSym);
    int i, n = length(Parent), pari, pos;
    int *ai, *ap, *nzc = Calloc(n, int), ntot;

    ntot = 0;
    for (i = n - 1; i >= 0; i--) { /* count non-zeros in each column */
    pari = INTEGER(Parent)[i];
    ntot += (nzc[i] = (pari >= 0) ? 1 + nzc[pari] : 0);
    }

    slot_dup(ans, x, Matrix_DimSym);
    slot_dup(ans, x, Matrix_DimNamesSym);
    slot_dup(ans, x, Matrix_uploSym);
    slot_dup(ans, x, Matrix_diagSym);
    ai = INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, ntot));
    ap = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
    ap[0] = pos = 0;
    for (i = 0; i < n; i++) {
    ap[i + 1] = ap[i] + nzc[i];
    for (pari = INTEGER(Parent)[i]; pari >= 0; pari = INTEGER(Parent)[pari])
        ai[pos++] = pari;
    }

    Free(nzc);
    UNPROTECT(1);
    return ans;
}

/** 
 * Solve  one  of the matrix equations  op(A)*C = B, or * C*op(A) = B
 * where A is an lCholCMatrix object and B and C are lgCMatrix 
 * objects.
 *
 * An early check is done to see if A is the identity, in which case
 * BIP is returned and bp is unmodified.
 * 
 * @param side LFT or RGT
 * @param transa TRN or NTR
 * @param m number of rows in B and C
 * @param n number of columns in B and C
 * @param Parent Parent array of A
 * @param BIP pointer to the row indices for B
 * @param bp array of column pointers for B (may be overwritten)
 *
 * @return Pointer to the updated row indices for C.  The column
 * pointers bp are overwritten with cp.
 */
SEXP
lCholClgCsm(enum CBLAS_SIDE side, enum CBLAS_TRANSPOSE transa, int m,
        int n, const int Parent[], SEXP BIP, int bp[])
{
    int *bi = INTEGER(BIP), bnz, extra, ident, j, pari, pos;
    int nca = (transa == TRN) ? n : m;

    ident = 1;
    for (j = 0; j < nca; j++)
    if (Parent[j] >= 0) {
        ident = 0;
        break;
    }
    if (ident) return BIP;

    extra = 0;
    if (side == LFT) {
    if (transa == TRN) {
        error(_("code not yet written"));
        return R_NilValue;
    } else {
        int *Tci, *Tj, *Ti, ntot;
        for (j = 0; j < n; j++) {
        int ii, ii2 = bp[j + 1];
        for (ii = bp[j]; ii < ii2; ii++)
            for (pari = Parent[bi[ii]]; pari >= 0;
             pari = Parent[pari]) extra++;
        }

        bnz = bp[n];
        ntot = bnz + extra;
        Ti = Memcpy(Calloc(ntot, int), bi, bnz);
        Tj = expand_cmprPt(n, bp, Calloc(ntot, int));
        Tci = Calloc(ntot, int);

        pos = bnz;
        for (j = 0; j < n; j++) {
        int ii, ii2 = bp[j + 1];
        for (ii = bp[j]; ii < ii2; ii++)
            for (pari = Parent[bi[ii]]; pari >= 0;
             pari = Parent[pari]) {
                Ti[pos] = pari;
                Tj[pos] = j;
                pos++;
            }
        }
        triplet_to_col(m, n, ntot, Ti, Tj, (double *) NULL,
               bp, Tci, (double *) NULL);

        bnz = bp[n];
        BIP = PROTECT(allocVector(INTSXP, bnz));
        Memcpy(INTEGER(BIP), Tci, bnz);
        
        Free(Tci); Free(Ti); Free(Tj);
        UNPROTECT(1);
        return BIP;
    }
    } else {
    if (transa == TRN) {
        int *Tci, *Tj, *Ti, ntot;
        for (j = 0; j < n; j++) {
        int ii, ii2 = bp[j + 1];
        for (pari = Parent[j]; pari >= 0; pari = Parent[pari])
            for (ii = bp[j]; ii < ii2; ii++) extra++;
        }

        bnz = bp[n];
        ntot = bnz + extra;
        Ti = Memcpy(Calloc(ntot, int), bi, bnz);
        Tj = expand_cmprPt(n, bp, Calloc(ntot, int));
        Tci = Calloc(ntot, int);

        pos = bnz;
        for (j = 0; j < n; j++) {
        int ii, ii2 = bp[j + 1];
        for (pari = Parent[j]; pari >= 0; pari = Parent[pari]) {
            for (ii = bp[j]; ii < ii2; ii++) {
            Ti[pos] = bi[ii];
            Tj[pos] = pari;
            pos++;
            }
        }
        }
        triplet_to_col(m, n, ntot, Ti, Tj, (double *) NULL,
               bp, Tci, (double *) NULL);
        bnz = bp[n];
        BIP = PROTECT(allocVector(INTSXP, bnz));
        Memcpy(INTEGER(BIP), Tci, bnz);
        
        Free(Tci); Free(Ti); Free(Tj);
        UNPROTECT(1);
        return BIP;
    } else {
        error(_("code not yet written"));
        return R_NilValue;
    }
    }
}

SEXP lCholCMatrix_lgCMatrix_solve(SEXP a, SEXP b)
{
    SEXP ans = PROTECT(duplicate(b));
    int n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0],
    *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym));

    if (n != bdims[0])
    error(_("Dimensions of a (%d,%d) and b (%d,%d) do not conform"),
          n, n, bdims[0], bdims[1]);
    SET_SLOT(ans, Matrix_iSym,
         lCholClgCsm(LFT, NTR, bdims[0], bdims[1],
             INTEGER(GET_SLOT(a, Matrix_ParentSym)),
             GET_SLOT(ans, Matrix_iSym),
             INTEGER(GET_SLOT(ans, Matrix_pSym))));
    UNPROTECT(1);
    return ans;
}