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 2731, Mon Oct 17 18:07:09 2011 UTC
 Line 1
                          /* Sparse matrices in compressed column-oriented form */
  #include "Csparse.h"
  #include "Tsparse.h"
  #include "chm_common.h"
-Line 136
+Line 137
              }
      }
      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 155
      /* 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 176
                                GET_SLOT(x, Matrix_DimNamesSym));
  }
+ // n.CMatrix --> [dli].CMatrix  (not going through CHM!)
+ SEXP nz_pattern_to_Csparse(SEXP x, SEXP res_kind)
+ {
+     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 = strdup(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)
  {
-     return chm_dense_to_matrix(cholmod_l_sparse_to_dense(AS_CHM_SP__(x), &c),
+     return chm_dense_to_matrix(cholmod_sparse_to_dense(AS_CHM_SP__(x), &c),
 /*do_free*/, GET_SLOT(x, Matrix_DimNamesSym));
  }
  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 258
      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 216
+Line 271
      int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
      R_CheckStack();
-     chgx = cholmod_l_copy(chx, /* stype: */ uploT, chx->xtype, &c);
+     chgx = cholmod_copy(chx, /* stype: */ uploT, chx->xtype, &c);
      /* xtype: pattern, "real", complex or .. */
      return chm_sparse_to_SEXP(chgx, 1, 0, Rkind, "",
                                GET_SLOT(x, Matrix_DimNamesSym));
-Line 228
+Line 283
       *       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 247
+Line 302
      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 353
      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 327
+Line 382
      CHM_SP cha = AS_CHM_SP(a);
      SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b));
      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);
  }
-Line 347
+Line 413
      CHM_SP cha = AS_CHM_SP(a);
      SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b));
      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 446
  #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 473
      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 503
      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 516
      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 524
  {
      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 542
      }
      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 529
+Line 601
      if (csize >= 0 && !isInteger(j))
          error(_("Index j must be NULL or integer"));
-     if (chx->stype) /* symmetricMatrix */
+     if (!chx->stype) {/* non-symmetric Matrix */
-         /* for now, cholmod_submatrix() only accepts "generalMatrix" */
+         return chm_sparse_to_SEXP(cholmod_submatrix(chx,
-         chx = cholmod_l_copy(chx, /* stype: */ 0, chx->xtype, &c);
-     return chm_sparse_to_SEXP(cholmod_l_submatrix(chx,
                                  (rsize < 0) ? NULL : INTEGER(i), rsize,
                                  (csize < 0) ? NULL : INTEGER(j), csize,
                                                    TRUE, TRUE, &c),
 , 0, Rkind, "",
                                /* FIXME: drops dimnames */ R_NilValue);
  }
+                                 /* for now, cholmod_submatrix() only accepts "generalMatrix" */
+     CHM_SP tmp = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);
+     CHM_SP ans = cholmod_submatrix(tmp,
+                                    (rsize < 0) ? NULL : INTEGER(i), rsize,
+                                    (csize < 0) ? NULL : INTEGER(j), csize,
+                                    TRUE, TRUE, &c);
+     cholmod_free_sparse(&tmp, &c);
+     return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", 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)
  {
-Line 548
+Line 643
      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 764
+Line 859
      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 869
          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 887
      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.2731
 Legend:



Removed from v.2586
 


changed lines


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

root@r-forge.r-project.org	ViewVC Help
Powered by ViewVC 1.0.0