SCM

SCM Repository

[matrix] View of /branches/Matrix-mer2/src/lCholCMatrix.c
ViewVC logotype

View of /branches/Matrix-mer2/src/lCholCMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 733 - (download) (as text) (annotate)
Tue May 17 17:18:30 2005 UTC (14 years, 9 months ago) by bates
Original Path: pkg/src/lCholCMatrix.c
File size: 5787 byte(s)
Adding logical sparse matrix methods
				/* 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.
 * 
 * @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 bi array of row indices of B
 * @param bp array of column pointers of B
 * @param CIP pointer to the row indices for C
 * @param cp array of column pointers for C
 *
 * @return Pointer to the updated row indices for C.  The column
 * pointers may also be overwritten.
 */
SEXP
lCholClgCsm(enum CBLAS_SIDE side, enum CBLAS_TRANSPOSE transa, int m,
	    int n, const int Parent[], const int bi[], const int bp[],
	    SEXP CIP, int cp[])
{
    int *ci = INTEGER(CIP), cnz, extra, i, j, pari, pos;

    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 = bi[ii]; pari >= 0; pari = Parent[pari])
			if (check_csc_index(cp, ci, pari, j, -1) < 0) extra++;
	    }
/* This is not quite right.  C may contain more elements than in the solution. */
	    if (!extra) return CIP;

	    cnz = cp[n];
	    ntot = cnz + extra;
	    Ti = Memcpy(Calloc(ntot, int), ci, cnz);
	    Tj = expand_cmprPt(n, cp, Calloc(ntot, int));
	    Tci = Calloc(ntot, int);

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

	    cnz = cp[n];
	    CIP = PROTECT(allocVector(INTSXP, cnz));
	    Memcpy(INTEGER(CIP), Tci, cnz);
	    
	    Free(Tci); Free(Ti); Free(Tj);
	    UNPROTECT(1);
	    return CIP;
	}
    } else {
	if (transa == TRN) {
	    int *Tci, *Tj, *Ti, ntot;
	    for (j = 0; j < n; j++) {
		int ii, ii2 = bp[j + 1];
		for (pari = j; pari >= 0; pari = Parent[pari]) {
		    for (ii = bp[j]; ii < ii2; ii++)
			if (check_csc_index(cp, ci, ii, pari, -1) < 0) extra++;
	    }
/* This is not quite right.  C may contain more elements than in the solution. */
	    if (!extra) return CIP;

	    cnz = cp[n];
	    ntot = cnz + extra;
	    Ti = Memcpy(Calloc(ntot, int), ci, cnz);
	    Tj = expand_cmprPt(n, cp, Calloc(ntot, int));
	    Tci = Calloc(ntot, int);

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

	    cnz = cp[n];
	    CIP = PROTECT(allocVector(INTSXP, cnz));
	    Memcpy(INTEGER(CIP), Tci, cnz);
	    
	    Free(Tci); Free(Ti); Free(Tj);
	    UNPROTECT(1);
	    return CIP;
	} 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)),
			 INTEGER(GET_SLOT(b, Matrix_iSym)),
			 INTEGER(GET_SLOT(b, Matrix_pSym)),
			 GET_SLOT(ans, Matrix_iSym),
			 INTEGER(GET_SLOT(ans, Matrix_pSym))));
    UNPROTECT(1);
    return ans;
}
			 
	     

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge