Logo Search packages:      
Sourcecode: rmatrix version File versions  Download package

lCholCMatrix.c

                        /* 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;
}
                   
           

Generated by  Doxygen 1.6.0   Back to index