SCM Repository

[matrix] Diff of /pkg/Matrix/src/Csparse.c

[matrix] / pkg / Matrix / src / Csparse.c

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

-revision 2586, Sun Jul 25 02:32:06 2010 UTC
+revision 3020, Tue Oct 14 16:14:02 2014 UTC
 Line 1
                          /* Sparse matrices in compressed column-oriented form */
  #include "Csparse.h"
  #include "Tsparse.h"
  #include "chm_common.h"
-Line 36
+Line 37
      return Csparse_validate_(x, FALSE);
  }
- SEXP Csparse_validate2(SEXP x, SEXP maybe_modify) {
-     return Csparse_validate_(x, asLogical(maybe_modify));
- }
- SEXP Csparse_validate_(SEXP x, Rboolean maybe_modify)
+ #define _t_Csparse_validate
- {
+ #include "t_Csparse_validate.c"
-     /* NB: we do *NOT* check a potential 'x' slot here, at all */
-     SEXP pslot = GET_SLOT(x, Matrix_pSym),
-         islot = GET_SLOT(x, Matrix_iSym);
-     Rboolean sorted, strictly;
-     int j, k,
-         *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
-         nrow = dims[0],
-         ncol = dims[1],
-         *xp = INTEGER(pslot),
-         *xi = INTEGER(islot);
-     if (length(pslot) != dims[1] + 1)
+ #define _t_Csparse_sort
-         return mkString(_("slot p must have length = ncol(.) + 1"));
+ #include "t_Csparse_validate.c"
-     if (xp[0] != 0)
-         return mkString(_("first element of slot p must be zero"));
-     if (length(islot) < xp[ncol]) /* allow larger slots from over-allocation!*/
-         return
-             mkString(_("last element of slot p must match length of slots i and x"));
-     for (j = 0; j < xp[ncol]; j++) {
-         if (xi[j] < 0 || xi[j] >= nrow)
-             return mkString(_("all row indices must be between 0 and nrow-1"));
-     }
-     sorted = TRUE; strictly = TRUE;
-     for (j = 0; j < ncol; j++) {
-         if (xp[j] > xp[j + 1])
-             return mkString(_("slot p must be non-decreasing"));
-         if(sorted) /* only act if >= 2 entries in column j : */
-             for (k = xp[j] + 1; k < xp[j + 1]; k++) {
-                 if (xi[k] < xi[k - 1])
-                     sorted = FALSE;
-                 else if (xi[k] == xi[k - 1])
-                     strictly = FALSE;
-             }
-     }
-     if (!sorted) {
-         if(maybe_modify) {
-             CHM_SP chx = (CHM_SP) alloca(sizeof(cholmod_sparse));
-             R_CheckStack();
-             as_cholmod_sparse(chx, x, FALSE, TRUE);/*-> cholmod_l_sort() ! */
-             /* as chx = AS_CHM_SP__(x)  but  ^^^^ sorting x in_place !!! */
-             /* Now re-check that row indices are *strictly* increasing
+ // R: .validateCsparse(x, sort.if.needed = FALSE) :
-              * (and not just increasing) within each column : */
+ SEXP Csparse_validate2(SEXP x, SEXP maybe_modify) {
-             for (j = 0; j < ncol; j++) {
+     return Csparse_validate_(x, asLogical(maybe_modify));
-                 for (k = xp[j] + 1; k < xp[j + 1]; k++)
-                     if (xi[k] == xi[k - 1])
-                         return mkString(_("slot i is not *strictly* increasing inside a column (even after cholmod_l_sort)"));
-             }
-         } else { /* no modifying sorting : */
-             return mkString(_("row indices are not sorted within columns"));
-         }
-     } else if(!strictly) {  /* sorted, but not strictly */
-         return mkString(_("slot i is not *strictly* increasing inside a column"));
      }
-     return ScalarLogical(1);
+ // R: Matrix:::.sortCsparse(x) :
+ SEXP Csparse_sort (SEXP x) {
+    int ok = Csparse_sort_2(x, TRUE); // modifying x directly
+    if(!ok) warning(_("Csparse_sort(x): x is not a valid (apart from sorting) CsparseMatrix"));
+    return x;
  }
  SEXP Rsparse_validate(SEXP x)
-Line 136
+Line 93
              }
      }
      if (!sorted)
-         /* cannot easily use cholmod_l_sort(.) ... -> "error out" :*/
+         /* cannot easily use cholmod_sort(.) ... -> "error out" :*/
          return mkString(_("slot j is not increasing inside a column"));
      else if(!strictly) /* sorted, but not strictly */
          return mkString(_("slot j is not *strictly* increasing inside a column"));
-Line 154
+Line 111
      /* This loses the symmetry property, since cholmod_dense has none,
       * BUT, much worse (FIXME!), it also transforms CHOLMOD_PATTERN ("n") matrices
       * to numeric (CHOLMOD_REAL) ones : */
-     CHM_DN chxd = cholmod_l_sparse_to_dense(chxs, &c);
+     CHM_DN chxd = cholmod_sparse_to_dense(chxs, &c);
      int Rkind = (chxs->xtype == CHOLMOD_PATTERN)? -1 : Real_kind(x);
      R_CheckStack();
      return chm_dense_to_SEXP(chxd, 1, Rkind, GET_SLOT(x, Matrix_DimNamesSym));
  }
+ // FIXME: do not go via CHM (should not be too hard, to just *drop* the x-slot, right?
  SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri)
  {
      CHM_SP chxs = AS_CHM_SP__(x);
-     CHM_SP chxcp = cholmod_l_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c);
+     CHM_SP chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c);
      int tr = asLogical(tri);
      R_CheckStack();
-Line 174
+Line 132
                                GET_SLOT(x, Matrix_DimNamesSym));
  }
- SEXP Csparse_to_matrix(SEXP x)
+ // n.CMatrix --> [dli].CMatrix  (not going through CHM!)
+ SEXP nz_pattern_to_Csparse(SEXP x, SEXP res_kind)
  {
-     return chm_dense_to_matrix(cholmod_l_sparse_to_dense(AS_CHM_SP__(x), &c),
+     return nz2Csparse(x, asInteger(res_kind));
+ }
+ // n.CMatrix --> [dli].CMatrix  (not going through CHM!)
+ SEXP nz2Csparse(SEXP x, enum x_slot_kind r_kind)
+ {
+     const char *cl_x = class_P(x);
+     if(cl_x[0] != 'n') error(_("not a 'n.CMatrix'"));
+     if(cl_x[2] != 'C') error(_("not a CsparseMatrix"));
+     int nnz = LENGTH(GET_SLOT(x, Matrix_iSym));
+     SEXP ans;
+     char *ncl = alloca(strlen(cl_x) + 1); /* not much memory required */
+     strcpy(ncl, cl_x);
+     double *dx_x; int *ix_x;
+     ncl[0] = (r_kind == x_double ? 'd' :
+               (r_kind == x_logical ? 'l' :
+                /* else (for now):  r_kind == x_integer : */ 'i'));
+     PROTECT(ans = NEW_OBJECT(MAKE_CLASS(ncl)));
+     // create a correct 'x' slot:
+     switch(r_kind) {
+         int i;
+     case x_double: // 'd'
+         dx_x = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz));
+         for (i=0; i < nnz; i++) dx_x[i] = 1.;
+         break;
+     case x_logical: // 'l'
+         ix_x = LOGICAL(ALLOC_SLOT(ans, Matrix_xSym, LGLSXP, nnz));
+         for (i=0; i < nnz; i++) ix_x[i] = TRUE;
+         break;
+     case x_integer: // 'i'
+         ix_x = INTEGER(ALLOC_SLOT(ans, Matrix_xSym, INTSXP, nnz));
+         for (i=0; i < nnz; i++) ix_x[i] = 1;
+         break;
+     default:
+         error(_("nz2Csparse(): invalid/non-implemented r_kind = %d"),
+               r_kind);
+     }
+     // now copy all other slots :
+     slot_dup(ans, x, Matrix_iSym);
+     slot_dup(ans, x, Matrix_pSym);
+     slot_dup(ans, x, Matrix_DimSym);
+     slot_dup(ans, x, Matrix_DimNamesSym);
+     if(ncl[1] != 'g') { // symmetric or triangular ...
+         slot_dup_if_has(ans, x, Matrix_uploSym);
+         slot_dup_if_has(ans, x, Matrix_diagSym);
+     }
+     UNPROTECT(1);
+     return ans;
+ }
+ SEXP Csparse_to_matrix(SEXP x, SEXP chk)
+ {
+     return chm_dense_to_matrix(cholmod_sparse_to_dense(AS_CHM_SP2(x, asLogical(chk)), &c),
 /*do_free*/, GET_SLOT(x, Matrix_DimNamesSym));
  }
+ SEXP Csparse_to_vector(SEXP x)
+ {
+     return chm_dense_to_vector(cholmod_sparse_to_dense(AS_CHM_SP__(x), &c), 1);
+ }
  SEXP Csparse_to_Tsparse(SEXP x, SEXP tri)
  {
      CHM_SP chxs = AS_CHM_SP__(x);
-     CHM_TR chxt = cholmod_l_sparse_to_triplet(chxs, &c);
+     CHM_TR chxt = cholmod_sparse_to_triplet(chxs, &c);
      int tr = asLogical(tri);
      int Rkind = (chxs->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
      R_CheckStack();
-Line 203
+Line 219
      if (!(chx->stype))
          error(_("Nonsymmetric matrix in Csparse_symmetric_to_general"));
-     chgx = cholmod_l_copy(chx, /* stype: */ 0, chx->xtype, &c);
+     chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);
      /* xtype: pattern, "real", complex or .. */
      return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",
                                GET_SLOT(x, Matrix_DimNamesSym));
-Line 211
+Line 227
  SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo)
  {
+     int *adims = INTEGER(GET_SLOT(x, Matrix_DimSym)), n = adims[0];
+     if(n != adims[1]) {
+         error(_("Csparse_general_to_symmetric(): matrix is not square!"));
+         return R_NilValue; /* -Wall */
+     }
      CHM_SP chx = AS_CHM_SP__(x), chgx;
-     int uploT = (*CHAR(STRING_ELT(uplo,0)) == 'U') ? 1 : -1;
+     int uploT = (*CHAR(asChar(uplo)) == 'U') ? 1 : -1;
      int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
      R_CheckStack();
+     chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c);
-     chgx = cholmod_l_copy(chx, /* stype: */ uploT, chx->xtype, &c);
+     /* need _symmetric_ dimnames */
+     SEXP dns = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))),
+         nms_dns = getAttrib(dns, R_NamesSymbol);
+     if(!equal_string_vectors(VECTOR_ELT(dns, 0),
+                              VECTOR_ELT(dns, 1))) {
+         if(uploT == 1)
+             SET_VECTOR_ELT(dns, 0, VECTOR_ELT(dns,1));
+         else
+             SET_VECTOR_ELT(dns, 1, VECTOR_ELT(dns,0));
+     }
+     if(!isNull(nms_dns) &&  // names(dimnames(.)) :
+        !R_compute_identical(STRING_ELT(nms_dns, 0),
+                             STRING_ELT(nms_dns, 1), 15)) {
+         if(uploT == 1)
+             SET_STRING_ELT(nms_dns, 0, STRING_ELT(nms_dns,1));
+         else
+             SET_STRING_ELT(nms_dns, 1, STRING_ELT(nms_dns,0));
+         setAttrib(dns, R_NamesSymbol, nms_dns);
+     }
+     UNPROTECT(1);
      /* xtype: pattern, "real", complex or .. */
-     return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",
+     return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "", dns);
-                               GET_SLOT(x, Matrix_DimNamesSym));
  }
  SEXP Csparse_transpose(SEXP x, SEXP tri)
-Line 228
+Line 269
       *       since cholmod (& cs) lacks sparse 'int' matrices */
      CHM_SP chx = AS_CHM_SP__(x);
      int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
-     CHM_SP chxt = cholmod_l_transpose(chx, chx->xtype, &c);
+     CHM_SP chxt = cholmod_transpose(chx, chx->xtype, &c);
      SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp;
      int tr = asLogical(tri);
      R_CheckStack();
-Line 236
+Line 277
      tmp = VECTOR_ELT(dn, 0);    /* swap the dimnames */
      SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1));
      SET_VECTOR_ELT(dn, 1, tmp);
+     if(!isNull(tmp = getAttrib(dn, R_NamesSymbol))) { // swap names(dimnames(.)):
+         SEXP nms_dns = PROTECT(allocVector(VECSXP, 2));
+         SET_VECTOR_ELT(nms_dns, 1, STRING_ELT(tmp, 0));
+         SET_VECTOR_ELT(nms_dns, 0, STRING_ELT(tmp, 1));
+         setAttrib(dn, R_NamesSymbol, nms_dns);
+         UNPROTECT(1);
+     }
      UNPROTECT(1);
      return chm_sparse_to_SEXP(chxt, 1, /* SWAP 'uplo' for triangular */
                                tr ? ((*uplo_P(x) == 'U') ? -1 : 1) : 0,
-Line 247
+Line 295
      CHM_SP
          cha = AS_CHM_SP(a),
          chb = AS_CHM_SP(b),
-         chc = cholmod_l_ssmult(cha, chb, /*out_stype:*/ 0,
+         chc = cholmod_ssmult(cha, chb, /*out_stype:*/ 0,
                                 /* values:= is_numeric (T/F) */ cha->xtype > 0,
                                 /*out sorted:*/ 1, &c);
      const char *cl_a = class_P(a), *cl_b = class_P(b);
-Line 298
+Line 346
      SEXP dn = PROTECT(allocVector(VECSXP, 2));
      R_CheckStack();
-     chTr = cholmod_l_transpose((tr) ? chb : cha, chb->xtype, &c);
+     chTr = cholmod_transpose((tr) ? chb : cha, chb->xtype, &c);
-     chc = cholmod_l_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb,
+     chc = cholmod_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb,
                           /*out_stype:*/ 0, cha->xtype, /*out sorted:*/ 1, &c);
-     cholmod_l_free_sparse(&chTr, &c);
+     cholmod_free_sparse(&chTr, &c);
      /* Preserve triangularity and unit-triangularity if appropriate;
       * see Csparse_Csparse_prod() for comments */
-Line 325
+Line 373
  SEXP Csparse_dense_prod(SEXP a, SEXP b)
  {
      CHM_SP cha = AS_CHM_SP(a);
-     SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b));
+     SEXP b_M = PROTECT(mMatrix_as_dgeMatrix2(b, // transpose_if_vector =
+                                              cha->ncol == 1));
      CHM_DN chb = AS_CHM_DN(b_M);
-     CHM_DN chc = cholmod_l_allocate_dense(cha->nrow, chb->ncol, cha->nrow,
+     CHM_DN chc = cholmod_allocate_dense(cha->nrow, chb->ncol, cha->nrow,
                                          chb->xtype, &c);
      SEXP dn = PROTECT(allocVector(VECSXP, 2));
      double one[] = {1,0}, zero[] = {0,0};
+     int nprot = 2;
      R_CheckStack();
+     /* Tim Davis, please FIXME:  currently (2010-11) *fails* when  a  is a pattern matrix:*/
-     cholmod_l_sdmult(cha, 0, one, zero, chb, chc, &c);
+     if(cha->xtype == CHOLMOD_PATTERN) {
+         /* warning(_("Csparse_dense_prod(): cholmod_sdmult() not yet implemented for pattern./ ngCMatrix" */
+         /*        " --> slightly inefficient coercion")); */
+         // This *fails* to produce a CHOLMOD_REAL ..
+         // CHM_SP chd = cholmod_l_copy(cha, cha->stype, CHOLMOD_REAL, &c);
+         // --> use our Matrix-classes
+         SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++;
+         cha = AS_CHM_SP(da);
+     }
+     cholmod_sdmult(cha, 0, one, zero, chb, chc, &c);
      SET_VECTOR_ELT(dn, 0,       /* establish dimnames */
                     duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0)));
      SET_VECTOR_ELT(dn, 1,
                     duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1)));
-     UNPROTECT(2);
+     UNPROTECT(nprot);
      return chm_dense_to_SEXP(chc, 1, 0, dn);
  }
  SEXP Csparse_dense_crossprod(SEXP a, SEXP b)
  {
      CHM_SP cha = AS_CHM_SP(a);
-     SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b));
+     SEXP b_M = PROTECT(mMatrix_as_dgeMatrix2(b, // transpose_if_vector =
+                                              cha->nrow == 1));
      CHM_DN chb = AS_CHM_DN(b_M);
-     CHM_DN chc = cholmod_l_allocate_dense(cha->ncol, chb->ncol, cha->ncol,
+     CHM_DN chc = cholmod_allocate_dense(cha->ncol, chb->ncol, cha->ncol,
                                          chb->xtype, &c);
-     SEXP dn = PROTECT(allocVector(VECSXP, 2));
+     SEXP dn = PROTECT(allocVector(VECSXP, 2)); int nprot = 2;
      double one[] = {1,0}, zero[] = {0,0};
      R_CheckStack();
+     // -- see Csparse_dense_prod() above :
-     cholmod_l_sdmult(cha, 1, one, zero, chb, chc, &c);
+     if(cha->xtype == CHOLMOD_PATTERN) {
+         SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++;
+         cha = AS_CHM_SP(da);
+     }
+     cholmod_sdmult(cha, 1, one, zero, chb, chc, &c);
      SET_VECTOR_ELT(dn, 0,       /* establish dimnames */
                     duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1)));
      SET_VECTOR_ELT(dn, 1,
                     duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1)));
-     UNPROTECT(2);
+     UNPROTECT(nprot);
      return chm_dense_to_SEXP(chc, 1, 0, dn);
  }
-Line 376
+Line 441
  #endif
      CHM_SP chcp, chxt,
          chx = (trip ?
-                cholmod_l_triplet_to_sparse(cht, cht->nnz, &c) :
+                cholmod_triplet_to_sparse(cht, cht->nnz, &c) :
                 AS_CHM_SP(x));
      SEXP dn = PROTECT(allocVector(VECSXP, 2));
      R_CheckStack();
-     if (!tr) chxt = cholmod_l_transpose(chx, chx->xtype, &c);
+     if (!tr) chxt = cholmod_transpose(chx, chx->xtype, &c);
-     chcp = cholmod_l_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c);
+     chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c);
      if(!chcp) {
          UNPROTECT(1);
-         error(_("Csparse_crossprod(): error return from cholmod_l_aat()"));
+         error(_("Csparse_crossprod(): error return from cholmod_aat()"));
      }
-     cholmod_l_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c);
+     cholmod_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c);
      chcp->stype = 1;
-     if (trip) cholmod_l_free_sparse(&chx, &c);
+     if (trip) cholmod_free_sparse(&chx, &c);
-     if (!tr) cholmod_l_free_sparse(&chxt, &c);
+     if (!tr) cholmod_free_sparse(&chxt, &c);
      SET_VECTOR_ELT(dn, 0,       /* establish dimnames */
                     duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym),
                                          (tr) ? 0 : 1)));
-Line 403
+Line 468
      return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn);
  }
+ /* Csparse_drop(x, tol):  drop entries with absolute value < tol, i.e,
+ *  at least all "explicit" zeros */
  SEXP Csparse_drop(SEXP x, SEXP tol)
  {
      const char *cl = class_P(x);
      /* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */
      int tr = (cl[1] == 't');
      CHM_SP chx = AS_CHM_SP__(x);
-     CHM_SP ans = cholmod_l_copy(chx, chx->stype, chx->xtype, &c);
+     CHM_SP ans = cholmod_copy(chx, chx->stype, chx->xtype, &c);
      double dtol = asReal(tol);
      int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
      R_CheckStack();
-     if(!cholmod_l_drop(dtol, ans, &c))
+     if(!cholmod_drop(dtol, ans, &c))
-         error(_("cholmod_l_drop() failed"));
+         error(_("cholmod_drop() failed"));
      return chm_sparse_to_SEXP(ans, 1,
                                tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
                                Rkind, tr ? diag_P(x) : "",
-Line 431
+Line 498
      R_CheckStack();
      /* TODO: currently drops dimnames - and we fix at R level */
-     return chm_sparse_to_SEXP(cholmod_l_horzcat(chx, chy, 1, &c),
+     return chm_sparse_to_SEXP(cholmod_horzcat(chx, chy, 1, &c),
 , 0, Rkind, "", R_NilValue);
  }
-Line 444
+Line 511
      R_CheckStack();
      /* TODO: currently drops dimnames - and we fix at R level */
-     return chm_sparse_to_SEXP(cholmod_l_vertcat(chx, chy, 1, &c),
+     return chm_sparse_to_SEXP(cholmod_vertcat(chx, chy, 1, &c),
 , 0, Rkind, "", R_NilValue);
  }
-Line 452
+Line 519
  {
      CHM_SP chx = AS_CHM_SP__(x);
      int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
-     CHM_SP ans = cholmod_l_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c);
+     CHM_SP ans = cholmod_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c);
      R_CheckStack();
      return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "",
-Line 470
+Line 537
      }
      else { /* unit triangular (diag='U'): "fill the diagonal" & diag:= "N" */
          CHM_SP chx = AS_CHM_SP__(x);
-         CHM_SP eye = cholmod_l_speye(chx->nrow, chx->ncol, chx->xtype, &c);
+         CHM_SP eye = cholmod_speye(chx->nrow, chx->ncol, chx->xtype, &c);
          double one[] = {1, 0};
-         CHM_SP ans = cholmod_l_add(chx, eye, one, one, TRUE, TRUE, &c);
+         CHM_SP ans = cholmod_add(chx, eye, one, one, TRUE, TRUE, &c);
          int uploT = (*uplo_P(x) == 'U') ? 1 : -1;
          int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
          R_CheckStack();
-         cholmod_l_free_sparse(&eye, &c);
+         cholmod_free_sparse(&eye, &c);
          return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N",
                                    GET_SLOT(x, Matrix_DimNamesSym));
      }
-Line 494
+Line 561
      }
      else { /* triangular with diag='N'): now drop the diagonal */
          /* duplicate, since chx will be modified: */
-         CHM_SP chx = AS_CHM_SP__(duplicate(x));
+         SEXP xx = PROTECT(duplicate(x));
+         CHM_SP chx = AS_CHM_SP__(xx);
          int uploT = (*uplo_P(x) == 'U') ? 1 : -1,
              Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
          R_CheckStack();
          chm_diagN2U(chx, uploT, /* do_realloc */ FALSE);
-         return chm_sparse_to_SEXP(chx, /*dofree*/ 0/* or 1 ?? */,
+         SEXP ans = chm_sparse_to_SEXP(chx, /*dofree*/ 0/* or 1 ?? */,
                                    uploT, Rkind, "U",
                                    GET_SLOT(x, Matrix_DimNamesSym));
+         UNPROTECT(1);// only now !
+         return ans;
      }
  }
-Line 529
+Line 599
      if (csize >= 0 && !isInteger(j))
          error(_("Index j must be NULL or integer"));
-     if (chx->stype) /* symmetricMatrix */
+ #define CHM_SUB(_M_, _i_, _j_)                                  \
+     cholmod_submatrix(_M_,                                      \
+                       (rsize < 0) ? NULL : INTEGER(_i_), rsize, \
+                       (csize < 0) ? NULL : INTEGER(_j_), csize, \
+                       TRUE, TRUE, &c)
+     CHM_SP ans;
+     if (!chx->stype) {/* non-symmetric Matrix */
+         ans = CHM_SUB(chx, i, j);
+     }
+     else {
          /* for now, cholmod_submatrix() only accepts "generalMatrix" */
-         chx = cholmod_l_copy(chx, /* stype: */ 0, chx->xtype, &c);
+         CHM_SP tmp = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);
+         ans = CHM_SUB(tmp, i, j);
+         cholmod_free_sparse(&tmp, &c);
+     }
-     return chm_sparse_to_SEXP(cholmod_l_submatrix(chx,
+     // "FIXME": currently dropping dimnames, and adding them afterwards in R :
-                                 (rsize < 0) ? NULL : INTEGER(i), rsize,
+     /* // dimnames: */
-                                 (csize < 0) ? NULL : INTEGER(j), csize,
+     /* SEXP x_dns = GET_SLOT(x, Matrix_DimNamesSym), */
-                                                   TRUE, TRUE, &c),
+     /*  dn = PROTECT(allocVector(VECSXP, 2)); */
-, 0, Rkind, "",
+     return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", /* dimnames: */ R_NilValue);
-                               /* FIXME: drops dimnames */ R_NilValue);
  }
+ #define _d_Csp_
+ #include "t_Csparse_subassign.c"
+ #define _l_Csp_
+ #include "t_Csparse_subassign.c"
+ #define _i_Csp_
+ #include "t_Csparse_subassign.c"
+ #define _n_Csp_
+ #include "t_Csparse_subassign.c"
+ #define _z_Csp_
+ #include "t_Csparse_subassign.c"
  SEXP Csparse_MatrixMarket(SEXP x, SEXP fname)
  {
      FILE *f = fopen(CHAR(asChar(fname)), "w");
-Line 548
+Line 646
      if (!f)
          error(_("failure to open file \"%s\" for writing"),
                CHAR(asChar(fname)));
-     if (!cholmod_l_write_sparse(f, AS_CHM_SP(x),
+     if (!cholmod_write_sparse(f, AS_CHM_SP(x),
                                (CHM_SP)NULL, (char*) NULL, &c))
-         error(_("cholmod_l_write_sparse returned error code"));
+         error(_("cholmod_write_sparse returned error code"));
      fclose(f);
      return R_NilValue;
  }
-Line 568
+Line 666
   *
   * @return  a SEXP, either a (double) number or a length n-vector of diagonal entries
   */
- SEXP diag_tC_ptr(int n, int *x_p, double *x_x, int *perm, SEXP resultKind)
+ SEXP diag_tC_ptr(int n, int *x_p, double *x_x, Rboolean is_U, int *perm,
  /*                                ^^^^^^ FIXME[Generalize] to int / ... */
+                  SEXP resultKind)
  {
      const char* res_ch = CHAR(STRING_ELT(resultKind,0));
-     enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log
+     enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log, min, max, range
      } res_kind = ((!strcmp(res_ch, "trace")) ? trace :
                    ((!strcmp(res_ch, "sumLog")) ? sum_log :
                     ((!strcmp(res_ch, "prod")) ? prod :
+                     ((!strcmp(res_ch, "min")) ? min :
+                      ((!strcmp(res_ch, "max")) ? max :
+                       ((!strcmp(res_ch, "range")) ? range :
                      ((!strcmp(res_ch, "diag")) ? diag :
                       ((!strcmp(res_ch, "diagBack")) ? diag_backpermuted :
-                       -1)))));
+                          -1))))))));
-     int i, n_x, i_from = 0;
+     int i, n_x, i_from;
      SEXP ans = PROTECT(allocVector(REALSXP,
  /*                                 ^^^^  FIXME[Generalize] */
                                     (res_kind == diag ||
-                                     res_kind == diag_backpermuted) ? n : 1));
+                                     res_kind == diag_backpermuted) ? n :
+                                    (res_kind == range ? 2 : 1)));
      double *v = REAL(ans);
  /*  ^^^^^^      ^^^^  FIXME[Generalize] */
+     i_from = (is_U ? -1 : 0);
  #define for_DIAG(v_ASSIGN)                                              \
-     for(i = 0; i < n; i++, i_from += n_x) {                             \
+     for(i = 0; i < n; i++) {                                    \
          /* looking at i-th column */                                    \
          n_x = x_p[i+1] - x_p[i];/* #{entries} in this column */ \
+         if( is_U) i_from += n_x;                                \
          v_ASSIGN;                                                       \
+         if(!is_U) i_from += n_x;                                \
      }
      /* NOTA BENE: we assume  -- uplo = "L" i.e. lower triangular matrix
       *            for uplo = "U" (makes sense with a "dtCMatrix" !),
-      *            should use  x_x[i_from + (nx - 1)] instead of x_x[i_from],
+      *            should use  x_x[i_from + (n_x - 1)] instead of x_x[i_from],
-      *            where nx = (x_p[i+1] - x_p[i])
+      *            where n_x = (x_p[i+1] - x_p[i])
       */
      switch(res_kind) {
-     case trace:
+     case trace: // = sum
          v[0] = 0.;
          for_DIAG(v[0] += x_x[i_from]);
          break;
-Line 616
+Line 723
          for_DIAG(v[0] *= x_x[i_from]);
          break;
+     case min:
+         v[0] = R_PosInf;
+         for_DIAG(if(v[0] > x_x[i_from]) v[0] = x_x[i_from]);
+         break;
+     case max:
+         v[0] = R_NegInf;
+         for_DIAG(if(v[0] < x_x[i_from]) v[0] = x_x[i_from]);
+         break;
+     case range:
+         v[0] = R_PosInf;
+         v[1] = R_NegInf;
+         for_DIAG(if(v[0] > x_x[i_from]) v[0] = x_x[i_from];
+                  if(v[1] < x_x[i_from]) v[1] = x_x[i_from]);
+         break;
      case diag:
          for_DIAG(v[i] = x_x[i_from]);
          break;
-Line 623
+Line 747
      case diag_backpermuted:
          for_DIAG(v[i] = x_x[i_from]);
-         warning(_("resultKind = 'diagBack' (back-permuted) is experimental"));
+         warning(_("%s = '%s' (back-permuted) is experimental"),
+                 "resultKind", "diagBack");
          /* now back_permute : */
          for(i = 0; i < n; i++) {
              double tmp = v[i]; v[i] = v[perm[i]]; v[perm[i]] = tmp;
-Line 644
+Line 769
   * Extract the diagonal entries from *triangular* Csparse matrix  __or__ a
   * cholmod_sparse factor (LDL = TRUE).
   *
+  * @param obj -- now a cholmod_sparse factor or a dtCMatrix
   * @param pslot  'p' (column pointer)   slot of Csparse matrix/factor
   * @param xslot  'x' (non-zero entries) slot of Csparse matrix/factor
   * @param perm_slot  'perm' (= permutation vector) slot of corresponding CHMfactor;
-Line 652
+Line 778
   *
   * @return  a SEXP, either a (double) number or a length n-vector of diagonal entries
   */
- SEXP diag_tC(SEXP pslot, SEXP xslot, SEXP perm_slot, SEXP resultKind)
+ SEXP diag_tC(SEXP obj, SEXP resultKind)
  {
+     SEXP
+         pslot = GET_SLOT(obj, Matrix_pSym),
+         xslot = GET_SLOT(obj, Matrix_xSym);
+     Rboolean is_U = (R_has_slot(obj, Matrix_uploSym) &&
+                      *CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))) == 'U');
      int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */
-         *x_p  = INTEGER(pslot),
+         *x_p  = INTEGER(pslot), pp = -1, *perm;
-         *perm = INTEGER(perm_slot);
      double *x_x = REAL(xslot);
  /*  ^^^^^^        ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/
-     return diag_tC_ptr(n, x_p, x_x, perm, resultKind);
+     if(R_has_slot(obj, Matrix_permSym))
+         perm = INTEGER(GET_SLOT(obj, Matrix_permSym));
+     else perm = &pp;
+     return diag_tC_ptr(n, x_p, x_x, is_U, perm, resultKind);
  }
  /**
   * Create a Csparse matrix object from indices and/or pointers.
   *
-Line 764
+Line 900
      if (cls[1] != 'g')
          error(_("Only 'g'eneral sparse matrix types allowed"));
                                  /* allocate and populate the triplet */
-     T = cholmod_l_allocate_triplet((size_t)nrow, (size_t)ncol, (size_t)nnz, 0,
+     T = cholmod_allocate_triplet((size_t)nrow, (size_t)ncol, (size_t)nnz, 0,
                                      xtype, &c);
      T->x = x;
      tri = (int*)T->i;
-Line 774
+Line 910
          trj[ii] = j[ii] - ((!mj && index1) ? 1 : 0);
      }
                                  /* create the cholmod_sparse structure */
-     A = cholmod_l_triplet_to_sparse(T, nnz, &c);
+     A = cholmod_triplet_to_sparse(T, nnz, &c);
-     cholmod_l_free_triplet(&T, &c);
+     cholmod_free_triplet(&T, &c);
                                  /* copy the information to the SEXP */
      ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cls)));
  /* FIXME: This has been copied from chm_sparse_to_SEXP in chm_common.c */
                                  /* allocate and copy common slots */
-     nnz = cholmod_l_nnz(A, &c);
+     nnz = cholmod_nnz(A, &c);
      dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
      dims[0] = A->nrow; dims[1] = A->ncol;
      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, A->ncol + 1)), (int*)A->p, A->ncol + 1);
-Line 792
+Line 928
      case 'l':
          error(_("code not yet written for cls = \"lgCMatrix\""));
      }
-     cholmod_l_free_sparse(&A, &c);
+ /* FIXME: dimnames are *NOT* put there yet (if non-NULL) */
+     cholmod_free_sparse(&A, &c);
      UNPROTECT(1);
      return ans;
  }

 Legend:



Removed from v.2586
 


changed lines


 
Added in v.3020
 Legend:



Removed from v.2586
 


changed lines


 
Added in v.3020
-Removed from v.2586
+Added in v.3020

R-Forge@R-project.org	ViewVC Help
Powered by ViewVC 1.0.0