SCM

SCM Repository

[matrix] View of /pkg/Matrix/src/Csparse.c
ViewVC logotype

View of /pkg/Matrix/src/Csparse.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 925 - (download) (as text) (annotate)
Mon Sep 19 19:01:31 2005 UTC (13 years, 7 months ago) by maechler
Original Path: pkg/src/Csparse.c
File size: 3208 byte(s)
too many changes: rbind2 for dM*, cbind2(), "Math", "Arith" for sparse ones; more  "s/gTMatrix/TsparseMatrix/" fixes; quite a few more tests, incl some for *all* classes
			/* Sparse matrices in compress column-oriented form */
#include "Csparse.h"
#ifdef USE_CHOLMOD
#include "chm_common.h"
#endif	/* USE_CHOLMOD */

SEXP Csparse_validate(SEXP x)
{
    SEXP pslot = GET_SLOT(x, Matrix_pSym),
	islot = GET_SLOT(x, Matrix_iSym);
    int j, ncol = length(pslot) - 1,
	*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
	nrow, *xp = INTEGER(pslot),
	*xi = INTEGER(islot);

    nrow = dims[0];
    if (length(pslot) <= 0)
	return mkString(_("slot p must have length > 0"));
    if (xp[0] != 0)
	return mkString(_("first element of slot p must be zero"));
    if (length(islot) != xp[ncol])
	return mkString(_("last element of slot p must match length of slots i and x"));
    for (j = 0; j < ncol; j++) {
	if (xp[j] > xp[j+1])
	    return mkString(_("slot p must be non-decreasing"));
    }
    for (j = 0; j < length(islot); j++) {
	if (xi[j] < 0 || xi[j] >= nrow)
	    return mkString(_("all row indices must be between 0 and nrow-1"));
    }
    return ScalarLogical(1);
}

SEXP Csparse_to_Tsparse(SEXP x)
{
#ifdef USE_CHOLMOD
    cholmod_sparse *chxs = as_cholmod_sparse(x);
    cholmod_triplet *chxt = cholmod_sparse_to_triplet(chxs, &c);

    Free(chxs);
    return chm_triplet_to_SEXP(chxt, 1);
#else
    error("General conversion requires CHOLMOD");
    return R_NilValue;		/* -Wall */
#endif	/* USE_CHOLMOD */
}

SEXP Csparse_transpose(SEXP x)
{
#ifdef USE_CHOLMOD
    cholmod_sparse *chx = as_cholmod_sparse(x);
    cholmod_sparse *chxt = cholmod_transpose(chx, (int) chx->xtype, &c);

    Free(chx);
    return chm_sparse_to_SEXP(chxt, 1);
#else
    error("General conversion requires CHOLMOD");
    return R_NilValue;		/* -Wall */
#endif	/* USE_CHOLMOD */
}


SEXP Csparse_Csparse_prod(SEXP a, SEXP b)
{
#ifdef USE_CHOLMOD
    cholmod_sparse *cha = as_cholmod_sparse(a), *chb = as_cholmod_sparse(b);
    cholmod_sparse *chc = cholmod_ssmult(cha, chb, 0, (int) cha->xtype, 1, &c);

    Free(cha); Free(chb);
    return chm_sparse_to_SEXP(chc, 1);
#else
    error("General multiplication requires CHOLMOD");
    return R_NilValue;		/* -Wall */
#endif	/* USE_CHOLMOD */
}

SEXP Csparse_dense_prod(SEXP a, SEXP b)
{
#ifdef USE_CHOLMOD
    cholmod_sparse *cha = as_cholmod_sparse(a);
    cholmod_dense *chb = as_cholmod_dense(b);
    cholmod_dense *chc = cholmod_allocate_dense(cha->nrow, chb->ncol,
						cha->nrow, chb->xtype, &c);
    double alpha = 1, beta = 0;

    cholmod_sdmult(cha, 0, &alpha, &beta, chb, chc, &c);
    Free(cha); Free(chb);
    return chm_dense_to_SEXP(chc, 1);
#else
    error("General multiplication requires CHOLMOD");
    return R_NilValue;		/* -Wall */
#endif	/* USE_CHOLMOD */
}

SEXP Csparse_crossprod(SEXP x, SEXP trans)
{
#ifdef USE_CHOLMOD
    int tr = asLogical(trans);	/* gets reversed because _aat is trcrossprod */
    cholmod_sparse *chx = as_cholmod_sparse(x), *chcp, *chxt;

    if (!tr)
	chxt = cholmod_transpose(chx, (int) chx->xtype, &c);
    chcp = cholmod_aat((tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c);

    Free(chx);
    if (!tr) cholmod_free_sparse(&chxt, &c);
    return chm_sparse_to_SEXP(chcp, 1);
#else
    error("General transpose requires CHOLMOD");
    return R_NilValue;		/* -Wall */
#endif	/* USE_CHOLMOD */
}


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