dropzero.c

Go to the documentation of this file.

#include "Mdefines.h"
#include "M5.h"
 
SEXP sparse_dropzero(SEXP from, const char *class, double tol)
{
    /* defined in ./aggregate.c : */
    SEXP sparse_aggregate(SEXP, const char *);
    from = sparse_aggregate(from, class);
    if (class[0] == 'n')
        return from;
    PROTECT(from);
 
    char ct = '\0';
    if (class[1] == 's' && class[0] == 'z')
        ct = TRANS(from);
 
    SEXP to;
    int strict = ISNAN(tol) || tol <= 0.0;
 
#define NZ(c, x) ((strict) ? c##NOT_ZERO(x) : c##NOT_ZERO_TOL(x, tol))
 
    if (class[2] != 'T') {
 
        SEXP iSym = (class[2] == 'C') ? Matrix_iSym : Matrix_jSym,
            p0 = PROTECT(GET_SLOT(from, Matrix_pSym)),
            i0 = PROTECT(GET_SLOT(from,        iSym)),
            x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
        int *pp0 = INTEGER(p0), *pi0 = INTEGER(i0),
            n = (int) (XLENGTH(p0) - 1), j, k, kend,
            nnz0 = pp0[n], nnz1 = 0;
        pp0++;
 
        if (ct != 'C') {
 
#define TEMPLATE(c) \
            do { \
                c##TYPE *px0 = c##PTR(x0); \
                for (k = 0; k < nnz0; ++k) \
                    if (NZ(c, px0[k])) \
                        ++nnz1; \
            } while (0)
 
            SWITCH4(class[0], TEMPLATE);
 
#undef TEMPLATE
 
        } else {
 
            Rcomplex *px0 = COMPLEX(x0);
            for (j = 0, k = 0; j < n; ++j) {
                kend = pp0[j];
                while (k < kend) {
                    if ((pi0[k] == j) ? NZ(d, px0[k].r) : NZ(z, px0[k]))
                        ++nnz1;
                    ++k;
                }
            }
 
        }
 
        if (nnz1 == nnz0) {
            UNPROTECT(4); /* x0, i0, p0, from */
            return from;
        }
 
        PROTECT(to = newObject(class));
 
        SEXP p1 = PROTECT(Rf_allocVector(INTSXP, XLENGTH(p0))),
            i1 = PROTECT(Rf_allocVector(INTSXP, nnz1)),
            x1 = PROTECT(Rf_allocVector(TYPEOF(x0), nnz1));
        int *pp1 = INTEGER(p1), *pi1 = INTEGER(i1);
        *(pp1++) = 0;
        SET_SLOT(to, Matrix_pSym, p1);
        SET_SLOT(to,        iSym, i1);
        SET_SLOT(to, Matrix_xSym, x1);
 
        if (ct != 'C') {
 
#define TEMPLATE(c) \
            do { \
                c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
                for (j = 0, k = 0; j < n; ++j) { \
                    pp1[j] = pp1[j - 1]; \
                    kend = pp0[j]; \
                    while (k < kend) { \
                        if (NZ(c, px0[k])) { \
                            pi1[pp1[j]  ] = pi0[k]; \
                            px1[pp1[j]++] = px0[k]; \
                        } \
                        ++k; \
                    } \
                } \
            } while (0)
 
            SWITCH4(class[0], TEMPLATE);
 
#undef TEMPLATE
 
        } else {
 
            Rcomplex *px0 = COMPLEX(x0), *px1 = COMPLEX(x1);
            for (j = 0, k = 0; j < n; ++j) {
                pp1[j] = pp1[j - 1];
                kend = pp0[j];
                while (k < kend) {
                    if ((pi0[k] == j) ? NZ(d, px0[k].r) : NZ(z, px0[k])) {
                        pi1[pp1[j]  ] = pi0[k];
                        px1[pp1[j]++] = px0[k];
                    }
                    ++k;
                }
            }
 
        }
 
        UNPROTECT(7); /* x1, i1, p1, to, x0, i0, p0 */
 
    } else {
 
        SEXP i0 = PROTECT(GET_SLOT(from, Matrix_iSym)),
            j0 = PROTECT(GET_SLOT(from, Matrix_jSym)),
            x0 = PROTECT(GET_SLOT(from, Matrix_xSym));
        int *pi0 = INTEGER(i0), *pj0 = INTEGER(j0);
        R_xlen_t k, nnz0 = XLENGTH(x0), nnz1 = 0;
 
        if (ct != 'C') {
 
#define TEMPLATE(c) \
            do { \
                c##TYPE *px0 = c##PTR(x0); \
                for (k = 0; k < nnz0; ++k) \
                    if (NZ(c, px0[k])) \
                        ++nnz1; \
            } while (0)
 
            SWITCH4(class[0], TEMPLATE);
 
#undef TEMPLATE
 
        } else {
 
            Rcomplex *px0 = COMPLEX(x0);
            for (k = 0; k < nnz0; ++k)
                if ((pi0[k] == pj0[k]) ? NZ(d, px0[k].r) : NZ(z, px0[k]))
                    ++nnz1;
 
        }
 
        if (nnz1 == nnz0) {
            UNPROTECT(4); /* x0, j0, i0, from */
            return from;
        }
 
        PROTECT(to = newObject(class));
 
        SEXP i1 = PROTECT(Rf_allocVector(INTSXP, nnz1)),
            j1 = PROTECT(Rf_allocVector(INTSXP, nnz1)),
            x1 = PROTECT(Rf_allocVector(TYPEOF(x0), nnz1));
        int *pi1 = INTEGER(i1), *pj1 = INTEGER(j1);
        SET_SLOT(to, Matrix_iSym, i1);
        SET_SLOT(to, Matrix_jSym, j1);
        SET_SLOT(to, Matrix_xSym, x1);
 
        if (ct != 'C') {
 
#define TEMPLATE(c) \
            do { \
                c##TYPE *px0 = c##PTR(x0), *px1 = c##PTR(x1); \
                for (k = 0; k < nnz0; ++k) \
                    if (NZ(c, px0[k])) { \
                        *(pi1++) = pi0[k]; \
                        *(pj1++) = pj0[k]; \
                        *(px1++) = px0[k]; \
                    } \
            } while (0)
 
            SWITCH4(class[0], TEMPLATE);
 
#undef TEMPLATE
 
        } else {
 
            Rcomplex *px0 = COMPLEX(x0), *px1 = COMPLEX(x1);
            for (k = 0; k < nnz0; ++k)
                if ((pi0[k] == pj0[k]) ? NZ(d, px0[k].r) : NZ(z, px0[k])) {
                    *(pi1++) = pi0[k];
                    *(pj1++) = pj0[k];
                    *(px1++) = px0[k];
                }
 
        }
 
        UNPROTECT(7); /* x1, j1, i1, to, x0, j0, i0 */
 
    }
 
    PROTECT(to);
    int *pdim = DIM(from), m = pdim[0], n = pdim[1];
    SET_DIM(to, m, n);
    SET_DIMNAMES(to, 0, DIMNAMES(from, 0));
    if (class[1] != 'g' && UPLO(from) != 'U')
        SET_UPLO(to);
    if (class[1] == 's' && class[0] == 'z' && ct != 'C')
        SET_TRANS(to);
    if (class[1] == 't' && DIAG(from) != 'N')
        SET_DIAG(to);
    UNPROTECT(2); /* to, from */
    return to;
}
 
SEXP R_sparse_dropzero(SEXP s_from, SEXP s_tol)
{
    const char *class = Matrix_class(s_from, valid_sparse, 6, __func__);
 
    double tol;
    if (TYPEOF(s_tol) != REALSXP || LENGTH(s_tol) < 1 ||
        ISNAN(tol = REAL(s_tol)[0]))
        Rf_error(_("'%s' is not a number"), "tol");
 
    return sparse_dropzero(s_from, class, tol);
}