SCM Repository

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

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

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

-revision 1371, Mon Aug 14 21:08:58 2006 UTC
+revision 1657, Wed Nov  1 16:29:53 2006 UTC
 Line 4
  SEXP Csparse_validate(SEXP x)
  {
+     /* 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);
-     int j, ncol = length(pslot) - 1,
+     int j, k, ncol, nrow, sorted,
          *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
-         nrow, *xp = INTEGER(pslot),
+         *xp = INTEGER(pslot),
          *xi = INTEGER(islot);
      nrow = dims[0];
-     if (length(pslot) <= 0)
+     ncol = dims[1];
-         return mkString(_("slot p must have length > 0"));
+     if (length(pslot) != dims[1] + 1)
+         return mkString(_("slot p must have length = ncol(.) + 1"));
      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"));
+         return
+             mkString(_("last element of slot p must match length of slots i and x"));
+     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"));
+     }
+     sorted = TRUE;
      for (j = 0; j < ncol; j++) {
          if (xp[j] > xp[j+1])
              return mkString(_("slot p must be non-decreasing"));
+         for (k = xp[j] + 1; k < xp[j + 1]; k++)
+             if (xi[k] < xi[k - 1]) sorted = FALSE;
      }
-     for (j = 0; j < length(islot); j++) {
+     if (!sorted) {
-         if (xi[j] < 0 || xi[j] >= nrow)
+         cholmod_sparse *chx = as_cholmod_sparse(x);
-             return mkString(_("all row indices must be between 0 and nrow-1"));
+         cholmod_sort(chx, &c);
+         Free(chx);
      }
      return ScalarLogical(1);
  }
-Line 35
+Line 46
      cholmod_dense *chxd = cholmod_sparse_to_dense(chxs, &c);
      Free(chxs);
-     return chm_dense_to_SEXP(chxd, 1);
+     return chm_dense_to_SEXP(chxd, 1, Real_kind(x));
  }
- SEXP Csparse_to_logical(SEXP x, SEXP tri)
+ SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri)
  {
      cholmod_sparse *chxs = as_cholmod_sparse(x);
      cholmod_sparse
-Line 51
+Line 62
              -1 : 1;
          diag = CHAR(asChar(GET_SLOT(x, Matrix_diagSym)));
      }
-     return chm_sparse_to_SEXP(chxcp, 1, uploT, diag,
+     return chm_sparse_to_SEXP(chxcp, 1, uploT, 0, diag,
                                  GET_SLOT(x, Matrix_DimNamesSym));
  }
-Line 69
+Line 80
  {
      cholmod_sparse *chxs = as_cholmod_sparse(x);
      cholmod_triplet *chxt = cholmod_sparse_to_triplet(chxs, &c);
-     int uploT = 0; char *diag = "";
+     int uploT = 0;
+     char *diag = "";
+     int Rkind = (chxs->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0;
      Free(chxs);
      if (asLogical(tri)) {       /* triangular sparse matrices */
-         uploT = (strcmp(CHAR(asChar(GET_SLOT(x, Matrix_uploSym))), "U")) ?
+         uploT = (*uplo_P(x) == 'U') ? -1 : 1;
-             -1 : 1;
+         diag = diag_P(x);
-         diag = CHAR(asChar(GET_SLOT(x, Matrix_diagSym)));
      }
-     return chm_triplet_to_SEXP(chxt, 1, uploT, diag,
+     return chm_triplet_to_SEXP(chxt, 1, uploT, Rkind, diag,
                                 GET_SLOT(x, Matrix_DimNamesSym));
  }
+ /* this used to be called  sCMatrix_to_gCMatrix(..)   [in ./dsCMatrix.c ]: */
  SEXP Csparse_symmetric_to_general(SEXP x)
  {
      cholmod_sparse *chx = as_cholmod_sparse(x), *chgx;
+     int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0;
      if (!(chx->stype))
-         error(_("Nonsymmetric matrix in Csparse_symmeteric_to_general"));
+         error(_("Nonsymmetric matrix in Csparse_symmetric_to_general"));
-     chgx = cholmod_copy(chx, 0, chx->xtype, &c);
+     chgx = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);
+     /* xtype: pattern, "real", complex or .. */
+     Free(chx);
+     return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",
+                               GET_SLOT(x, Matrix_DimNamesSym));
+ }
+ SEXP Csparse_general_to_symmetric(SEXP x, SEXP uplo)
+ {
+     cholmod_sparse *chx = as_cholmod_sparse(x), *chgx;
+     int uploT = (*CHAR(asChar(uplo)) == 'U') ? -1 : 1;
+     int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0;
+     chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c);
+     /* xtype: pattern, "real", complex or .. */
      Free(chx);
-     return chm_sparse_to_SEXP(chgx, 1, 0, "",
+     return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",
                                GET_SLOT(x, Matrix_DimNamesSym));
  }
  SEXP Csparse_transpose(SEXP x, SEXP tri)
  {
      cholmod_sparse *chx = as_cholmod_sparse(x);
+     int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0;
      cholmod_sparse *chxt = cholmod_transpose(chx, (int) chx->xtype, &c);
      SEXP dn = PROTECT(duplicate(GET_SLOT(x, Matrix_DimNamesSym))), tmp;
      int uploT = 0; char *diag = "";
-Line 106
+Line 135
      SET_VECTOR_ELT(dn, 1, tmp);
      UNPROTECT(1);
      if (asLogical(tri)) {       /* triangular sparse matrices */
-         uploT = (strcmp(CHAR(asChar(GET_SLOT(x, Matrix_uploSym))), "U")) ?
+         uploT = (*uplo_P(x) == 'U') ? -1 : 1;
-: -1;             /* switch upper and lower for transpose */
+         diag = diag_P(x);
-         diag = CHAR(asChar(GET_SLOT(x, Matrix_diagSym)));
      }
-     return chm_sparse_to_SEXP(chxt, 1, uploT, diag, dn);
+     return chm_sparse_to_SEXP(chxt, 1, uploT, Rkind, diag, dn);
  }
  SEXP Csparse_Csparse_prod(SEXP a, SEXP b)
-Line 125
+Line 153
                     duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0)));
      SET_VECTOR_ELT(dn, 1,
                     duplicate(VECTOR_ELT(GET_SLOT(b, Matrix_DimNamesSym), 1)));
-     return chm_sparse_to_SEXP(chc, 1, 0, "", dn);
+     return chm_sparse_to_SEXP(chc, 1, 0, 0, "", dn);
+ }
+ SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b)
+ {
+     cholmod_sparse *cha = as_cholmod_sparse(a),
+         *chb = as_cholmod_sparse(b);
+     cholmod_sparse *chta = cholmod_transpose(cha, 1, &c);
+     cholmod_sparse *chc = cholmod_ssmult(chta, chb, 0, cha->xtype, 1, &c);
+     SEXP dn = allocVector(VECSXP, 2);
+     Free(cha); Free(chb); cholmod_free_sparse(&chta, &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, Matrix_DimNamesSym), 1)));
+     return chm_sparse_to_SEXP(chc, 1, 0, 0, "", dn);
  }
  SEXP Csparse_dense_prod(SEXP a, SEXP b)
  {
      cholmod_sparse *cha = as_cholmod_sparse(a);
-     cholmod_dense *chb = as_cholmod_dense(b);
+     cholmod_dense *chb = as_cholmod_dense(PROTECT(mMatrix_as_dgeMatrix(b)));
-     cholmod_dense *chc = cholmod_allocate_dense(cha->nrow, chb->ncol,
+     cholmod_dense *chc =
-                                                 cha->nrow, chb->xtype, &c);
+         cholmod_allocate_dense(cha->nrow, chb->ncol, cha->nrow, chb->xtype, &c);
-     double alpha = 1, beta = 0;
+     double alpha[] = {1,0}, beta[] = {0,0};
-     cholmod_sdmult(cha, 0, &alpha, &beta, chb, chc, &c);
+     cholmod_sdmult(cha, 0, alpha, beta, chb, chc, &c);
      Free(cha); Free(chb);
-     return chm_dense_to_SEXP(chc, 1);
+     UNPROTECT(1);
+     return chm_dense_to_SEXP(chc, 1, 0);
  }
  SEXP Csparse_dense_crossprod(SEXP a, SEXP b)
  {
      cholmod_sparse *cha = as_cholmod_sparse(a);
-     cholmod_dense *chb = as_cholmod_dense(b);
+     cholmod_dense *chb = as_cholmod_dense(PROTECT(mMatrix_as_dgeMatrix(b)));
-     cholmod_dense *chc = cholmod_allocate_dense(cha->ncol, chb->ncol,
+     cholmod_dense *chc =
-                                                 cha->ncol, chb->xtype, &c);
+         cholmod_allocate_dense(cha->ncol, chb->ncol, cha->ncol, chb->xtype, &c);
-     double alpha = 1, beta = 0;
+     double alpha[] = {1,0}, beta[] = {0,0};
-     cholmod_sdmult(cha, 1, &alpha, &beta, chb, chc, &c);
+     cholmod_sdmult(cha, 1, alpha, beta, chb, chc, &c);
      Free(cha); Free(chb);
-     return chm_dense_to_SEXP(chc, 1);
+     UNPROTECT(1);
+     return chm_dense_to_SEXP(chc, 1, 0);
  }
  SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet)
-Line 169
+Line 215
          chxt = cholmod_transpose(chx, chx->xtype, &c);
      chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0, chx->xtype, &c);
      if(!chcp)
-         error("Csparse_crossprod(): error return from cholmod_aat()");
+         error(_("Csparse_crossprod(): error return from cholmod_aat()"));
      cholmod_band_inplace(0, chcp->ncol, chcp->xtype, chcp, &c);
      chcp->stype = 1;
      if (trip) {
-Line 185
+Line 231
                                          (tr) ? 1 : 0)));
      SET_VECTOR_ELT(dn, 1, duplicate(VECTOR_ELT(dn, 0)));
      UNPROTECT(1);
-     return chm_sparse_to_SEXP(chcp, 1, 0, "", dn);
+     return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn);
+ }
+ SEXP Csparse_drop(SEXP x, SEXP tol)
+ {
+     cholmod_sparse *chx = as_cholmod_sparse(x),
+         *ans = cholmod_copy(chx, chx->stype, chx->xtype, &c);
+     double dtol = asReal(tol);
+     int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0;
+     if(!cholmod_drop(dtol, ans, &c))
+         error(_("cholmod_drop() failed"));
+     Free(chx);
+     /* FIXME: currently drops dimnames */
+     return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue);
  }
  SEXP Csparse_horzcat(SEXP x, SEXP y)
  {
      cholmod_sparse *chx = as_cholmod_sparse(x),
          *chy = as_cholmod_sparse(y), *ans;
+     int Rkind = 0; /* only for "d" - FIXME */
      ans = cholmod_horzcat(chx, chy, 1, &c);
      Free(chx); Free(chy);
      /* FIXME: currently drops dimnames */
-     return chm_sparse_to_SEXP(ans, 1, 0, "", R_NilValue);
+     return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue);
  }
  SEXP Csparse_vertcat(SEXP x, SEXP y)
  {
      cholmod_sparse *chx = as_cholmod_sparse(x),
          *chy = as_cholmod_sparse(y), *ans;
+     int Rkind = 0; /* only for "d" - FIXME */
      ans = cholmod_vertcat(chx, chy, 1, &c);
      Free(chx); Free(chy);
      /* FIXME: currently drops dimnames */
-     return chm_sparse_to_SEXP(ans, 1, 0, "", R_NilValue);
+     return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue);
  }
  SEXP Csparse_band(SEXP x, SEXP k1, SEXP k2)
  {
      cholmod_sparse *chx = as_cholmod_sparse(x), *ans;
+     int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0;
      ans = cholmod_band(chx, asInteger(k1), asInteger(k2), chx->xtype, &c);
      Free(chx);
-     return chm_sparse_to_SEXP(ans, 1, 0, "", R_NilValue);
+     return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue);
  }
  SEXP Csparse_diagU2N(SEXP x)
-Line 227
+Line 291
      cholmod_sparse *ans = cholmod_add(chx, eye, one, one, TRUE, TRUE, &c);
      int uploT = (strcmp(CHAR(asChar(GET_SLOT(x, Matrix_uploSym))), "U")) ?
          -1 : 1;
+     int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0;
      Free(chx); cholmod_free_sparse(&eye, &c);
-     return chm_sparse_to_SEXP(ans, 1, uploT, "N",
+     return chm_sparse_to_SEXP(ans, 1, uploT, Rkind, "N",
                                duplicate(GET_SLOT(x, Matrix_DimNamesSym)));
  }
-Line 238
+Line 303
      cholmod_sparse *chx = as_cholmod_sparse(x);
      int rsize = (isNull(i)) ? -1 : LENGTH(i),
          csize = (isNull(j)) ? -1 : LENGTH(j);
+     int Rkind = (chx->xtype == CHOLMOD_REAL) ? Real_kind(x) : 0;
      if (rsize >= 0 && !isInteger(i))
          error(_("Index i must be NULL or integer"));
-Line 246
+Line 312
      return chm_sparse_to_SEXP(cholmod_submatrix(chx, INTEGER(i), rsize,
                                                  INTEGER(j), csize,
                                                  TRUE, TRUE, &c),
-, 0, "", R_NilValue);
+, 0, Rkind, "", R_NilValue);
  }

 Legend:



Removed from v.1371
 


changed lines


 
Added in v.1657
 Legend:



Removed from v.1371
 


changed lines


 
Added in v.1657
-Removed from v.1371
+Added in v.1657

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