# SCM Repository

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

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

revision 3072, Fri Mar 27 15:10:48 2015 UTC revision 3147, Thu Oct 29 16:56:10 2015 UTC
# Line 1  Line 1
1                          /* Sparse matrices in compressed column-oriented form */  /** @file Csparse.c
2     * The "CsparseMatrix" class from R package Matrix:
3     *
4     * Sparse matrices in compressed column-oriented form
5     */
6  #include "Csparse.h"  #include "Csparse.h"
7  #include "Tsparse.h"  #include "Tsparse.h"
8  #include "chm_common.h"  #include "chm_common.h"
# Line 101  Line 104
104      return ScalarLogical(1);      return ScalarLogical(1);
105  }  }
106
107  /**  /** @brief From a CsparseMatrix, produce a dense one.
108   * From a CsparseMatrix, produce a dense one.   *
109   * Directly deals with symmetric, triangular and general.   * Directly deals with symmetric, triangular and general.
110   * Called from ../R/Csparse.R's  C2dense()   * Called from ../R/Csparse.R's  C2dense()
111   *   *
# Line 122  Line 125
125          ctype = 0; // <- default = "dgC"          ctype = 0; // <- default = "dgC"
126      static const char *valid[] = { MATRIX_VALID_Csparse, ""};      static const char *valid[] = { MATRIX_VALID_Csparse, ""};
127      if(is_sym_or_tri == NA_INTEGER) { // find if  is(x, "symmetricMatrix") :      if(is_sym_or_tri == NA_INTEGER) { // find if  is(x, "symmetricMatrix") :
128          ctype = Matrix_check_class_etc(x, valid);          ctype = R_check_class_etc(x, valid);
129          is_sym = (ctype % 3 == 1);          is_sym = (ctype % 3 == 1);
130          is_tri = (ctype % 3 == 2);          is_tri = (ctype % 3 == 2);
131      } else {      } else {
# Line 130  Line 133
133          is_tri = is_sym_or_tri < 0;          is_tri = is_sym_or_tri < 0;
134          // => both are FALSE  iff  is_.. == 0          // => both are FALSE  iff  is_.. == 0
135          if(is_sym || is_tri)          if(is_sym || is_tri)
136              ctype = Matrix_check_class_etc(x, valid);              ctype = R_check_class_etc(x, valid);
137      }      }
138      CHM_SP chxs = AS_CHM_SP__(x);// -> chxs->stype = +- 1 <==> symmetric      CHM_SP chxs = AS_CHM_SP__(x);// -> chxs->stype = +- 1 <==> symmetric
139      R_CheckStack();      R_CheckStack();
# Line 219  Line 222
222  SEXP nz2Csparse(SEXP x, enum x_slot_kind r_kind)  SEXP nz2Csparse(SEXP x, enum x_slot_kind r_kind)
223  {  {
224      const char *cl_x = class_P(x);      const char *cl_x = class_P(x);
225      if(cl_x[0] != 'n') error(_("not a 'n.CMatrix'"));      // quick check - if ok, fast
226      if(cl_x[2] != 'C') error(_("not a CsparseMatrix"));      if(cl_x[0] != 'n' || cl_x[2] != 'C') {
227            // e.g. class = "A", from  setClass("A", contains = "ngCMatrix")
228            static const char *valid[] = { MATRIX_VALID_nCsparse, ""};
229            int ctype = R_check_class_etc(x, valid);
230            if(ctype < 0)
231                error(_("not a 'n.CMatrix'"));
232            else // fine : get a valid  cl_x  class_P()-like string :
233                cl_x = valid[ctype];
234        }
235      int nnz = LENGTH(GET_SLOT(x, Matrix_iSym));      int nnz = LENGTH(GET_SLOT(x, Matrix_iSym));
236      SEXP ans;      SEXP ans;
237      char *ncl = alloca(strlen(cl_x) + 1); /* not much memory required */      char *ncl = alloca(strlen(cl_x) + 1); /* not much memory required */
# Line 269  Line 280
280      int is_sym = asLogical(symm);      int is_sym = asLogical(symm);
281      if(is_sym == NA_LOGICAL) { // find if  is(x, "symmetricMatrix") :      if(is_sym == NA_LOGICAL) { // find if  is(x, "symmetricMatrix") :
282          static const char *valid[] = { MATRIX_VALID_Csparse, ""};          static const char *valid[] = { MATRIX_VALID_Csparse, ""};
283          int ctype = Matrix_check_class_etc(x, valid);          int ctype = R_check_class_etc(x, valid);
284          is_sym = (ctype % 3 == 1);          is_sym = (ctype % 3 == 1);
285      }      }
286      return chm_dense_to_matrix(      return chm_dense_to_matrix(
# Line 409  Line 420
420                                Rkind, tr ? diag_P(x) : "", dn);                                Rkind, tr ? diag_P(x) : "", dn);
421  }  }
422
423  /* NOTA BENE:  cholmod_ssmult(A,B, ...) ->  ./CHOLMOD/MatrixOps/cholmod_ssmult.c  /** @brief  A %*% B  - for matrices of class CsparseMatrix (R package "Matrix")
424     *
425     * @param a
426     * @param b
427     * @param bool_arith
428     *
429     * @return
430     *
431     * NOTA BENE:  cholmod_ssmult(A,B, ...) ->  ./CHOLMOD/MatrixOps/cholmod_ssmult.c
432   * ---------  computes a patter*n* matrix __always_ when   * ---------  computes a patter*n* matrix __always_ when
433   * *one* of A or B is pattern*n*, because of this (line 73-74):   * *one* of A or B is pattern*n*, because of this (line 73-74):
434     ---------------------------------------------------------------------------     ---------------------------------------------------------------------------
# Line 418  Line 437
437     ---------------------------------------------------------------------------     ---------------------------------------------------------------------------
438   * ==> Often need to copy the patter*n* to a *l*ogical matrix first !!!   * ==> Often need to copy the patter*n* to a *l*ogical matrix first !!!
439   */   */

440  SEXP Csparse_Csparse_prod(SEXP a, SEXP b, SEXP bool_arith)  SEXP Csparse_Csparse_prod(SEXP a, SEXP b, SEXP bool_arith)
441  {  {
442      CHM_SP      CHM_SP
443          cha = AS_CHM_SP(a),          cha = AS_CHM_SP(a),
444          chb = AS_CHM_SP(b), chc;          chb = AS_CHM_SP(b), chc;
445      R_CheckStack();      R_CheckStack();
// const char *cl_a = class_P(a), *cl_b = class_P(b);
446      static const char *valid_tri[] = { MATRIX_VALID_tri_Csparse, "" };      static const char *valid_tri[] = { MATRIX_VALID_tri_Csparse, "" };
447      char diag[] = {'\0', '\0'};      char diag[] = {'\0', '\0'};
448      int uploT = 0, nprot = 1,      int uploT = 0, nprot = 1,
# Line 463  Line 480
480       * Note that in that case, the multiplication itself should happen       * Note that in that case, the multiplication itself should happen
481       * faster.  But there's no support for that in CHOLMOD */       * faster.  But there's no support for that in CHOLMOD */
482
483      if(Matrix_check_class_etc(a, valid_tri) >= 0 &&      if(R_check_class_etc(a, valid_tri) >= 0 &&
484         Matrix_check_class_etc(b, valid_tri) >= 0)         R_check_class_etc(b, valid_tri) >= 0)
485          if(*uplo_P(a) == *uplo_P(b)) { /* both upper, or both lower tri. */          if(*uplo_P(a) == *uplo_P(b)) { /* both upper, or both lower tri. */
486              uploT = (*uplo_P(a) == 'U') ? 1 : -1;              uploT = (*uplo_P(a) == 'U') ? 1 : -1;
487              if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */              if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */
# Line 484  Line 501
501      return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn);      return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn);
502  }  }
503
504  /* trans = FALSE:  crossprod(a,b)  /** @brief [t]crossprod (<Csparse>, <Csparse>)
505   * trans = TRUE : tcrossprod(a,b) */   *
506     * @param a a "CsparseMatrix" object
507     * @param b a "CsparseMatrix" object
508     * @param trans trans = FALSE:  crossprod(a,b)
509     *              trans = TRUE : tcrossprod(a,b)
510     * @param bool_arith logical (TRUE / NA / FALSE): Should boolean arithmetic be used.
511     *
512     * @return a CsparseMatrix, the (t)cross product of a and b.
513     */
514  SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b, SEXP trans, SEXP bool_arith)  SEXP Csparse_Csparse_crossprod(SEXP a, SEXP b, SEXP trans, SEXP bool_arith)
515  {  {
516      int tr = asLogical(trans), nprot = 1,      int tr = asLogical(trans), nprot = 1,
# Line 522  Line 547
547          if(!a_is_n && !b_is_n) {          if(!a_is_n && !b_is_n) {
548              // coerce 'a' to pattern              // coerce 'a' to pattern
549              SEXP da = PROTECT(Csparse2nz(a, /* tri = */              SEXP da = PROTECT(Csparse2nz(a, /* tri = */
550                                           Matrix_check_class_etc(a, valid_tri) >= 0)); nprot++;                                           R_check_class_etc(a, valid_tri) >= 0)); nprot++;
551              cha = AS_CHM_SP(da);              cha = AS_CHM_SP(da);
552              R_CheckStack();              R_CheckStack();
553              // a_is_n = TRUE;              // a_is_n = TRUE;
# Line 536  Line 561
561
562      /* Preserve triangularity and unit-triangularity if appropriate;      /* Preserve triangularity and unit-triangularity if appropriate;
563       * see Csparse_Csparse_prod() for comments */       * see Csparse_Csparse_prod() for comments */
564      if(Matrix_check_class_etc(a, valid_tri) >= 0 &&      if(R_check_class_etc(a, valid_tri) >= 0 &&
565         Matrix_check_class_etc(b, valid_tri) >= 0)         R_check_class_etc(b, valid_tri) >= 0)
566          if(*uplo_P(a) != *uplo_P(b)) { /* one 'U', the other 'L' */          if(*uplo_P(a) != *uplo_P(b)) { /* one 'U', the other 'L' */
567              uploT = (*uplo_P(b) == 'U') ? 1 : -1;              uploT = (*uplo_P(b) == 'U') ? 1 : -1;
568              if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */              if(*diag_P(a) == 'U' && *diag_P(b) == 'U') { /* return UNIT-triag. */
# Line 561  Line 586
586  /**  /**
587   * All (dense * sparse)  Matrix products and cross products   * All (dense * sparse)  Matrix products and cross products
588   *   *
589   *   f( f(<Csparse>)  %*%  f(<dense>) )   for  f in {t(), identity()}   *   f( f(<Csparse>)  %*%  f(<dense>) )   where  f ()  is either t () [tranpose] or the identity.
590   *   *
591   * @param a CsparseMatrix  (n x m)   * @param a CsparseMatrix  (n x m)
592   * @param b numeric vector, matrix, or denseMatrix (m x k) or (k x m)  if `transp` is '2' or 'B'   * @param b numeric vector, matrix, or denseMatrix (m x k) or (k x m)  if `transp` is '2' or 'B'
# Line 605  Line 630
630      /* repeating a "cheap part" of  mMatrix_as_dgeMatrix2(b, .)  to see if      /* repeating a "cheap part" of  mMatrix_as_dgeMatrix2(b, .)  to see if
631       * we have a vector that we might 'transpose_if_vector' : */       * we have a vector that we might 'transpose_if_vector' : */
632      static const char *valid[] = {"_NOT_A_CLASS_", MATRIX_VALID_ddense, ""};      static const char *valid[] = {"_NOT_A_CLASS_", MATRIX_VALID_ddense, ""};
633      /* int ctype = Matrix_check_class_etc(b, valid);      /* int ctype = R_check_class_etc(b, valid);
634       * if (ctype > 0)   /.* a ddenseMatrix object */       * if (ctype > 0)   /.* a ddenseMatrix object */
635      if (Matrix_check_class_etc(b, valid) < 0) {      if (R_check_class_etc(b, valid) < 0) {
636          // not a ddenseM*:  is.matrix() or vector:          // not a ddenseM*:  is.matrix() or vector:
637          b_is_vector = !isMatrix(b);          b_is_vector = !isMatrix(b);
638      }      }
# Line 686  Line 711
711  }  }
712
713
714  /* Computes   x'x  or  x x' -- *also* for Tsparse (triplet = TRUE)  /** @brief Computes   x'x  or  x x' -- *also* for Tsparse (triplet = TRUE)
715     see Csparse_Csparse_crossprod above for  x'y and x y' */      see Csparse_Csparse_crossprod above for  x'y and x y'
716    */
717  SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet, SEXP bool_arith)  SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet, SEXP bool_arith)
718  {  {
719      int tripl = asLogical(triplet),      int tripl = asLogical(triplet),
# Line 699  Line 725
725      SEXP xx = PROTECT(Tsparse_diagU2N(x));      SEXP xx = PROTECT(Tsparse_diagU2N(x));
726      CHM_TR cht = tripl ? AS_CHM_TR__(xx) : (CHM_TR) NULL; int nprot = 2;      CHM_TR cht = tripl ? AS_CHM_TR__(xx) : (CHM_TR) NULL; int nprot = 2;
727  #endif  #endif
728      CHM_SP chcp, chxt,      CHM_SP chcp, chxt, chxc,
729          chx = (tripl ?          chx = (tripl ?
730                 cholmod_triplet_to_sparse(cht, cht->nnz, &c) :                 cholmod_triplet_to_sparse(cht, cht->nnz, &c) :
731                 AS_CHM_SP(x));                 AS_CHM_SP(x));
# Line 707  Line 733
733      R_CheckStack();      R_CheckStack();
734      Rboolean      Rboolean
735          x_is_n = (chx->xtype == CHOLMOD_PATTERN),          x_is_n = (chx->xtype == CHOLMOD_PATTERN),
736            x_is_sym = chx->stype != 0,
737          force_num = (do_bool == FALSE);          force_num = (do_bool == FALSE);
738
739      if(x_is_n && force_num) {      if(x_is_n && force_num) {
# Line 719  Line 746
746          // coerce 'x' to pattern          // coerce 'x' to pattern
747          static const char *valid_tri[] = { MATRIX_VALID_tri_Csparse, "" };          static const char *valid_tri[] = { MATRIX_VALID_tri_Csparse, "" };
748          SEXP dx = PROTECT(Csparse2nz(x, /* tri = */          SEXP dx = PROTECT(Csparse2nz(x, /* tri = */
749                                       Matrix_check_class_etc(x, valid_tri) >= 0)); nprot++;                                       R_check_class_etc(x, valid_tri) >= 0)); nprot++;
750          chx = AS_CHM_SP(dx);          chx = AS_CHM_SP(dx);
751          R_CheckStack();          R_CheckStack();
752      }      }
753
754      if (!tr) chxt = cholmod_transpose(chx, chx->xtype, &c);      if (!tr) chxt = cholmod_transpose(chx, chx->xtype, &c);
755
756      chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0,      if (x_is_sym) // cholmod_aat() does not like symmetric
757                         /* mode: */ chx->xtype, &c);          chxc = cholmod_copy(tr ? chx : chxt, /* stype: */ 0,
758                                chx->xtype, &c);
759        // CHOLMOD/Core/cholmod_aat.c :
760        chcp = cholmod_aat(x_is_sym ? chxc : (tr ? chx : chxt),
761                           (int *) NULL, 0, /* mode: */ chx->xtype, &c);
762      if(!chcp) {      if(!chcp) {
763          UNPROTECT(1);          UNPROTECT(1);
764          error(_("Csparse_crossprod(): error return from cholmod_aat()"));          error(_("Csparse_crossprod(): error return from cholmod_aat()"));
# Line 745  Line 776
776      return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn);      return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn);
777  }  }
778
779  /* Csparse_drop(x, tol):  drop entries with absolute value < tol, i.e,  /** @brief Csparse_drop(x, tol):  drop entries with absolute value < tol, i.e,
780  *  at least all "explicit" zeros */   *  at least all "explicit" zeros. */
781  SEXP Csparse_drop(SEXP x, SEXP tol)  SEXP Csparse_drop(SEXP x, SEXP tol)
782  {  {
783      const char *cl = class_P(x);      const char *cl = class_P(x);
784      /* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */      /* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */
785      int tr = (cl[1] == 't');      int tr = (cl[1] == 't'); // FIXME - rather  R_check_class_etc(..)
786      CHM_SP chx = AS_CHM_SP__(x);      CHM_SP chx = AS_CHM_SP__(x);
787      CHM_SP ans = cholmod_copy(chx, chx->stype, chx->xtype, &c);      CHM_SP ans = cholmod_copy(chx, chx->stype, chx->xtype, &c);
788      double dtol = asReal(tol);      double dtol = asReal(tol);
# Line 766  Line 797
797                                GET_SLOT(x, Matrix_DimNamesSym));                                GET_SLOT(x, Matrix_DimNamesSym));
798  }  }
799
800    /** @brief Horizontal Concatenation -  cbind( <Csparse>,  <Csparse>)
801     */
802  SEXP Csparse_horzcat(SEXP x, SEXP y)  SEXP Csparse_horzcat(SEXP x, SEXP y)
803  {  {
804  #define CSPARSE_CAT(_KIND_)                                             \  #define CSPARSE_CAT(_KIND_)                                             \
# Line 798  Line 831
831                                1, 0, Rkind, "", R_NilValue);                                1, 0, Rkind, "", R_NilValue);
832  }  }
833
834    /** @brief Vertical Concatenation -  rbind( <Csparse>,  <Csparse>)
835     */
836  SEXP Csparse_vertcat(SEXP x, SEXP y)  SEXP Csparse_vertcat(SEXP x, SEXP y)
837  {  {
838      CSPARSE_CAT("vertcat");      CSPARSE_CAT("vertcat");
# Line 870  Line 905
905  }  }
906
907  /**  /**
908   * "Indexing" aka subsetting : Compute  x[i,j], also for vectors i and j   * Indexing aka subsetting : Compute  x[i,j], also for vectors i and j
909   * Working via CHOLMOD_submatrix, see ./CHOLMOD/MatrixOps/cholmod_submatrix.c   * Working via CHOLMOD_submatrix, see ./CHOLMOD/MatrixOps/cholmod_submatrix.c
910   * @param x CsparseMatrix   * @param x CsparseMatrix
911   * @param i row     indices (0-origin), or NULL (R's)   * @param i row     indices (0-origin), or NULL (R, not C)
912   * @param j columns indices (0-origin), or NULL   * @param j columns indices (0-origin), or NULL
913   *   *
914   * @return x[i,j]  still CsparseMatrix --- currently, this loses dimnames   * @return x[i,j]  still CsparseMatrix --- currently, this loses dimnames

Legend:
 Removed from v.3072 changed lines Added in v.3147