t_subassign.c

Go to the documentation of this file.

/*------ Definition of a template for [diln]Csparse_subassign(...) : *
 *                       --------     ~~~~~~~~~~~~~~~~~~~~~~
 * i.e., included several times from ./Csparse.c
 *                                   ~~~~~~~~~~~
 * _slot_kind : use the integer codes matching x_slot_kind in ./Csparse.c
 */
 
#ifdef _d_Csp_
 
# define Csparse_subassign dCsparse_subassign
# define x_CLASSES "dgCMatrix",/* 0 */ "dtCMatrix" /* 1 */
# define sparseVECTOR "dsparseVector"
# define slot_kind_x 0
# define _DOUBLE_x
# define _has_x_slot_
 
# undef _d_Csp_
 
#elif defined (_l_Csp_)
 
# define Csparse_subassign lCsparse_subassign
# define x_CLASSES "lgCMatrix",/* 0 */ "ltCMatrix" /* 1 */
# define sparseVECTOR "lsparseVector"
# define slot_kind_x 1
# define _LGL_x
# define _has_x_slot_
 
# undef _l_Csp_
 
#elif defined (_i_Csp_)
 
# define Csparse_subassign iCsparse_subassign
# define x_CLASSES "igCMatrix",/* 0 */ "itCMatrix" /* 1 */
# define sparseVECTOR "isparseVector"
# define slot_kind_x 2
# define _INT_x
# define _has_x_slot_
 
# undef _i_Csp_
 
#elif defined (_n_Csp_)
 
# define Csparse_subassign nCsparse_subassign
# define x_CLASSES "ngCMatrix",/* 0 */ "ntCMatrix" /* 1 */
# define sparseVECTOR "nsparseVector"
# define slot_kind_x -1
# define _INT_x
/* withOUT 'x' slot -- CARE! we assume that this is the *ONLY* case w/o x slot */
 
# undef _n_Csp_
 
#elif defined (_z_Csp_)
 
// #  error "zgC* not yet implemented"
# define Csparse_subassign zCsparse_subassign
# define x_CLASSES "zgCMatrix",/* 0 */ "ztCMatrix" /* 1 */
# define sparseVECTOR "zsparseVector"
# define slot_kind_x 3
# define _CPLX_x
# define _has_x_slot_
 
# undef _z_Csp_
 
#else
 
#  error "no valid  _[dilnz]gC_ option"
 
#endif
// -------------------------------------------------
 
#ifdef _DOUBLE_x
 
# define Type_x double
# define Type_x_0_init(_VAR_) double _VAR_ = 0.
# define Type_x_1_init(_VAR_) double _VAR_ = 1.
# define STYP_x REAL
# define SXP_x  REALSXP
#undef _DOUBLE_x
 
#elif defined (_LGL_x)
 
# define Type_x int
# define Type_x_0_init(_VAR_) int _VAR_ = 0
# define Type_x_1_init(_VAR_) int _VAR_ = 1
# define STYP_x LOGICAL
# define SXP_x  LGLSXP
#undef _LGL_x
 
#elif defined (_INT_x)
 
# define Type_x int
# define Type_x_0_init(_VAR_) int _VAR_ = 0
# define Type_x_1_init(_VAR_) int _VAR_ = 1
# define STYP_x INTEGER
# define SXP_x  INTSXP
#undef _INT_x
 
#elif defined (_CPLX_x)
 
# define Type_x Rcomplex
# define Type_x_0_init(_VAR_) Rcomplex _VAR_;  _VAR_.r = _VAR_.i = 0.
# define Type_x_1_init(_VAR_) Rcomplex _VAR_;  _VAR_.r = 1.; _VAR_.i = 0.
# define STYP_x COMPLEX
# define SXP_x  CPLXSXP
 
#else
#  error "invalid macro logic"
#endif
 
 
 
SEXP Csparse_subassign(SEXP x, SEXP i_, SEXP j_, SEXP value)
{
    static const char
    *valid_cM [] = { // the only ones, for "the moment". FIXME: extend (!)
    x_CLASSES,
    ""},
    // value: assume a  "dsparseVector" for now -- slots: (i, length, x)
    *valid_spv[] = { sparseVECTOR, // = "the one with the same slot-class"
             // all others: ctype_v   slot_kind
             "nsparseVector",// 1     -1
             "lsparseVector",// 2      1
             "isparseVector",// 3      2
             "dsparseVector",// 4      0
             "zsparseVector",// 5      3
             ""};
 
    int ctype_x = R_check_class_etc(x,     valid_cM),
    ctype_v = R_check_class_etc(value, valid_spv);
    if (ctype_x < 0)
    Rf_error(_("invalid class of 'x' in Csparse_subassign()"));
    if (ctype_v < 0)
    Rf_error(_("invalid class of 'value' in Csparse_subassign()"));
    Rboolean value_is_nsp = ctype_v == 1;
#ifndef _has_x_slot_ // i.e. "n.CMatrix" : sparseVECTOR == "nsparseVector"
    if(!value_is_nsp) value_is_nsp = (ctype_v == 0);
#endif
 
    SEXP
    islot   = GET_SLOT(x, Matrix_iSym),
    dimslot = GET_SLOT(x, Matrix_DimSym),
    i_cp = PROTECT(Rf_coerceVector(i_, INTSXP)),
    j_cp = PROTECT(Rf_coerceVector(j_, INTSXP));
    // for d.CMatrix and l.CMatrix  but not n.CMatrix:
 
    int *dims = INTEGER(dimslot),
    ncol = dims[1], /* nrow = dims[0], */
    *i = INTEGER(i_cp), len_i = LENGTH(i_cp),
    *j = INTEGER(j_cp), len_j = LENGTH(j_cp),
    k,
    nnz_x = LENGTH(islot);
    int nnz = nnz_x;
 
#define MATRIX_SUBASSIGN_VERBOSE
// Temporary hack for debugging --- remove eventually -- FIXME
#ifdef MATRIX_SUBASSIGN_VERBOSE
    Rboolean verbose = i[0] < 0;
    if(verbose) {
    i[0] = -i[0];
    REprintf("Csparse_subassign() x[i,j] <- val; x is \"%s\"; value \"%s\" is_nsp=%d\n",
         valid_cM[ctype_x], valid_spv[ctype_v], (int)value_is_nsp);
    }
#endif
 
    SEXP val_i_slot, val_x_slot;
    val_i_slot = PROTECT(Rf_coerceVector(GET_SLOT(value, Matrix_iSym), REALSXP));
    double *val_i = REAL(val_i_slot);
    int nnz_val = LENGTH(GET_SLOT(value, Matrix_iSym)), n_prot = 4;
    Type_x *val_x = NULL;
    if(!value_is_nsp) {
    if(ctype_v) { // matrix 'x' and 'value' are of different kinds
        switch((enum x_slot_kind) slot_kind_x) {
        case x_pattern:// "n"
        case x_logical:// "l"
        if(ctype_v >= 3)
            Rf_warning(_("x[] <- val: val is coerced to logical for \"%s\" x"),
                   valid_cM[ctype_x]);
        break;
        case x_integer:
        if(ctype_v >= 4)
            Rf_error(_("x[] <- val: val should be integer or logical, is coerced to integer, for \"%s\" x"),
                 valid_cM[ctype_x]);
        break;
        case x_double:
        case x_complex: // coercion should be tried (and fail for complex -> double) below
        break;
        default:
        Rf_error(_("programming error in Csparse_subassign() should never happen"));
        }
        // otherwise: "coerce" :  as(., <sparseVector>) :
        val_x_slot = PROTECT(Rf_coerceVector(GET_SLOT(value, Matrix_xSym), SXP_x)); n_prot++;
        val_x = STYP_x(val_x_slot);
    } else {
        val_x = STYP_x(           GET_SLOT(value, Matrix_xSym));
    }
    }
    int64_t len_val = (int64_t) Rf_asReal(GET_SLOT(value, Matrix_lengthSym));
    /* llen_i = (int64_t) len_i; */
 
    SEXP ans;
    /* Instead of simple "duplicate" ans = Rf_duplicate(x), build up: */
    // Assuming that ans will have the same basic Matrix type as x :
    ans = PROTECT(newObject(valid_cM[ctype_x]));
    SET_SLOT(ans, Matrix_DimSym,      Rf_duplicate(dimslot));
    SET_SLOT(ans, Matrix_DimNamesSym, Rf_duplicate(GET_SLOT(x, Matrix_DimNamesSym)));
    SET_SLOT(ans, Matrix_pSym,        Rf_duplicate(GET_SLOT(x, Matrix_pSym)));
    SEXP r_pslot = GET_SLOT(ans, Matrix_pSym);
    // and assign the i- and x- slots at the end, as they are potentially modified
    // not just in content, but also in their *length*
    int *rp = INTEGER(r_pslot),
    *ri = R_Calloc(nnz_x, int);       // to contain the final i - slot
    Memcpy(ri, INTEGER(islot), nnz_x);
    Type_x_0_init(z_ans);
    Type_x_1_init(one_ans);
#ifdef _has_x_slot_
    Type_x *rx = R_Calloc(nnz_x, Type_x); // to contain the final x - slot
    Memcpy(rx, STYP_x(GET_SLOT(x, Matrix_xSym)), nnz_x);
#endif
    // NB:  nnz_x : will always be the "current allocated length" of (i, x) slots
    // --   nnz   : the current *used* length; always   nnz <= nnz_x
 
    int jj, j_val = 0; // in "running" conceptionally through all value[i+ jj*len_i]
    // values, we are "below"/"before" the (j_val)-th non-zero one.
    // e.g. if value = (0,0,...,0), have nnz_val == 0, j_val must remain == 0
    int64_t ii_val;// == "running" index (i + jj*len_i) % len_val for value[]
    for(jj = 0, ii_val=0; jj < len_j; jj++) {
    int j__ = j[jj];
    /* int64_t j_l = jj * llen_i; */
    R_CheckUserInterrupt();
    for(int ii = 0; ii < len_i; ii++, ii_val++) {
        int i__ = i[ii], p1, p2;
        if(nnz_val && ii_val >= len_val) { // "recycle" indexing into value[]
        ii_val -= len_val; // = (ii + jj*len_i) % len_val
        j_val = 0;
        }
        int64_t ii_v1;//= ii_val + 1;
        Type_x v, /* := value[(ii + j_l) % len_val]
             = .sparseVector_sub((ii + j_l) % len_val,
             nnz_val, val_i, val_x, len_val)
              */
        M_ij;
        int ind;
        Rboolean have_entry = FALSE;
 
        // note that rp[]'s may have *changed* even when 'j' remained!
        // "FIXME": do this only *when* rp[] has changed
        p1 = rp[j__], p2 = rp[j__ + 1];
 
        // v :=  value[(ii + j_l) % len_val] = value[ii_val]
        v = z_ans;
        if(j_val < nnz_val) { // maybe find v := non-zero value[ii_val]
        ii_v1 = ii_val + 1;
        if(ii_v1 < val_i[j_val]) { // typical case: are still in zero-stretch
            // v = z_ans (== 0)
        } else if(ii_v1 == val_i[j_val]) { // have a match
            v = (value_is_nsp) ? one_ans : val_x[j_val];
            j_val++;// from now on, look at the next non-zero entry
        } else { //  ii_v1 > val_i[j_val]
            REprintf("programming thinko in Csparse_subassign(*, i=%d,j=%d): ii_v=%lld, v@i[j_val=%d]=%g\n",
                 i__,j__, (long long)ii_v1, j_val, val_i[j_val]);
            j_val++;// from now on, look at the next non-zero entry
        }
        }
        // --------------- M_ij := getM(i., j.) --------------------------------
        M_ij = z_ans; // as in  ./t_sparseVector.c
        for(ind = p1; ind < p2; ind++) {
        if(ri[ind] >= i__) {
            if(ri[ind] == i__) {
#ifdef _has_x_slot_
            M_ij = rx[ind];
#else
            M_ij = 1;
#endif
#ifdef MATRIX_SUBASSIGN_VERBOSE
            if(verbose)
                REprintf("have entry x[%d, %d] = %g\n", i__, j__,
#  ifdef _CPLX_x
                     (double)M_ij.r);
#  else
                         (double)M_ij);
#  endif
#endif
                have_entry = TRUE;
            } else { // ri[ind] > i__
#ifdef MATRIX_SUBASSIGN_VERBOSE
                if(verbose)
                REprintf("@i > i__ = %d --> ind-- = %d\n", i__, ind);
#endif
        }
        break;
        }
    }
 
    //-- R:  if(getM(i., j.) != (v <- getV(ii, jj)))
 
    // if(contents differ) ==> value needs to be changed :
#ifdef _CPLX_x
    if(M_ij.r != v.r || M_ij.i != v.i) {
#else
        if(M_ij != v) {
#endif
 
#ifdef MATRIX_SUBASSIGN_VERBOSE
        if(verbose)
            REprintf("setting x[%d, %d] <- %g", i__,j__,
#  ifdef _CPLX_x
                 (double) v.r);
#  else
                     (double) v);
#  endif
#endif
        // (otherwise: nothing to do):
        // setM(i__, j__, v)
        // ----------------------------------------------------------
 
#ifndef _has_x_slot_
        if(v == z_ans) {
        // Case I ----- remove x[i, j] = M_ij  which we know is *non*-zero
//  BUT it is more efficient (memory-management!) *NOT* to remove,
//  currently using  drop0() in R code
        // we know : have_entry = TRUE ;
        //  ri[ind] == i__; M_ij = rx[ind];
#ifdef MATRIX_SUBASSIGN_VERBOSE
        if(verbose)
            REprintf(" rm ind=%d\n", ind);
#endif
        // remove the 'ind'-th element from x@i and x@x :
        nnz-- ;
        for(k=ind; k < nnz; k++) {
            ri[k] = ri[k+1];
#ifdef _has_x_slot_
            rx[k] = rx[k+1];
#endif
        }
        for(k=j__ + 1; k <= ncol; k++) {
            rp[k] = rp[k] - 1;
        }
        }
        else
#endif
        if(have_entry) {
            // Case II ----- replace (non-empty) x[i,j] by v -------
#ifdef MATRIX_SUBASSIGN_VERBOSE
            if(verbose)
            REprintf(" repl.  ind=%d\n", ind);
#endif
#ifdef _has_x_slot_
            rx[ind] = v;
#endif
        } else {
            // Case III ----- v != 0 : insert v into "empty" x[i,j] ----
 
            // extend the  i  and  x  slot by one entry : ---------------------
            if(nnz+1 > nnz_x) { // need to reallocate:
#ifdef MATRIX_SUBASSIGN_VERBOSE
            if(verbose) REprintf(" R_Realloc()ing: nnz_x=%d", nnz_x);
#endif
            // do it "only" 1x,..4x at the very most increasing by the
            // nnz-length of "value":
            nnz_x += (1 + nnz_val / 4);
#ifdef MATRIX_SUBASSIGN_VERBOSE
            if(verbose) REprintf("(nnz_v=%d) --> %d ", nnz_val, nnz_x);
#endif
            // C doc on realloc() says that the old content is *preserve*d
            ri = R_Realloc(ri, (size_t) nnz_x, int);
#ifdef _has_x_slot_
            rx = R_Realloc(rx, (size_t) nnz_x, Type_x);
#endif
            }
            // 3) fill them ...
 
            int i1 = ind;
#ifdef MATRIX_SUBASSIGN_VERBOSE
            if(verbose)
            REprintf(" INSERT p12=(%d,%d) -> ind=%d -> i1 = %d\n",
                 p1,p2, ind, i1);
#endif
 
            // shift the "upper values" *before* the insertion:
            for(int l = nnz-1; l >= i1; l--) {
            ri[l+1] = ri[l];
#ifdef _has_x_slot_
            rx[l+1] = rx[l];
#endif
            }
            ri[i1] = i__;
#ifdef _has_x_slot_
            rx[i1] = v;
#endif
            nnz++;
            // the columns j "right" of the current one :
            for(k=j__ + 1; k <= ncol; k++)
            rp[k]++;
        }
    }
#ifdef MATRIX_SUBASSIGN_VERBOSE
    else if(verbose) REprintf("M_ij == v = %g\n",
#  ifdef _CPLX_x
                  (double) v.r);
#  else
                              (double) v);
#  endif
#endif
    }// for( ii )
    }// for( jj )
 
    if(ctype_x == 1) { // triangularMatrix: copy the 'diag' and 'uplo' slots
    SET_SLOT(ans, Matrix_uploSym, Rf_duplicate(GET_SLOT(x, Matrix_uploSym)));
    SET_SLOT(ans, Matrix_diagSym, Rf_duplicate(GET_SLOT(x, Matrix_diagSym)));
    }
    // now assign the i- and x- slots,  free memory and return :
    PROTECT(islot = Rf_allocVector(INTSXP, nnz));
    Memcpy(INTEGER(islot), ri, nnz);
    SET_SLOT(ans, Matrix_iSym, islot);
    UNPROTECT(1);
#ifdef _has_x_slot_
    PROTECT(islot = Rf_allocVector(SXP_x, nnz));
    Memcpy(STYP_x(islot), rx, nnz);
    SET_SLOT(ans, Matrix_xSym, islot);
    UNPROTECT(1);
    R_Free(rx);
#endif
    R_Free(ri);
    UNPROTECT(n_prot);
    return ans;
}
 
#undef Csparse_subassign
#undef x_CLASSES
#undef sparseVECTOR
 
#undef Type_x
#undef STYP_x
#undef SXP_x
#undef Type_x_0_init
#undef Type_x_1_init
#undef _has_x_slot_
#undef slot_kind_x
 
#ifdef _CPLX_x
# undef _CPLX_x
#endif