SCM

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3069, Thu Mar 26 10:00:49 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 186  Line 189 
189  }  }
190    
191  // FIXME: do not go via CHM (should not be too hard, to just *drop* the x-slot, right?  // FIXME: do not go via CHM (should not be too hard, to just *drop* the x-slot, right?
192  SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri)  SEXP Csparse2nz(SEXP x, Rboolean tri)
193  {  {
194      CHM_SP chxs = AS_CHM_SP__(x);      CHM_SP chxs = AS_CHM_SP__(x);
195      CHM_SP chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c);      CHM_SP chxcp = cholmod_copy(chxs, chxs->stype, CHOLMOD_PATTERN, &c);
     int tr = asLogical(tri);  
196      R_CheckStack();      R_CheckStack();
197    
198      return chm_sparse_to_SEXP(chxcp, 1/*do_free*/,      return chm_sparse_to_SEXP(chxcp, 1/*do_free*/,
199                                tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,                                tri ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
200                                /* Rkind: pattern */ 0,                                /* Rkind: pattern */ 0,
201                                /* diag = */ tr ? diag_P(x) : "",                                /* diag = */ tri ? diag_P(x) : "",
202                                GET_SLOT(x, Matrix_DimNamesSym));                                GET_SLOT(x, Matrix_DimNamesSym));
203  }  }
204    SEXP Csparse_to_nz_pattern(SEXP x, SEXP tri)
205    {
206        int tr_ = asLogical(tri);
207        if(tr_ == NA_LOGICAL) {
208            warning(_("Csparse_to_nz_pattern(x, tri = NA): 'tri' is taken as TRUE"));
209            tr_ = TRUE;
210        }
211        return Csparse2nz(x, (Rboolean) tr_);
212    }
213    
214  // n.CMatrix --> [dli].CMatrix  (not going through CHM!)  // n.CMatrix --> [dli].CMatrix  (not going through CHM!)
215  SEXP nz_pattern_to_Csparse(SEXP x, SEXP res_kind)  SEXP nz_pattern_to_Csparse(SEXP x, SEXP res_kind)
# Line 211  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 261  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 401  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 410  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 455  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 476  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 487  Line 520 
520          chb = AS_CHM_SP(b),          chb = AS_CHM_SP(b),
521          chTr, chc;          chTr, chc;
522      R_CheckStack();      R_CheckStack();
     // const char *cl_a = class_P(a), *cl_b = class_P(b);  
523      static const char *valid_tri[] = { MATRIX_VALID_tri_Csparse, "" };      static const char *valid_tri[] = { MATRIX_VALID_tri_Csparse, "" };
524      char diag[] = {'\0', '\0'};      char diag[] = {'\0', '\0'};
525      int uploT = 0;      int uploT = 0;
# Line 502  Line 534 
534          SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++;          SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++;
535          cha = AS_CHM_SP(da);          cha = AS_CHM_SP(da);
536          R_CheckStack();          R_CheckStack();
537          a_is_n = FALSE;          // a_is_n = FALSE;
538      }      }
539      else if(b_is_n && (force_num || (maybe_bool && !a_is_n))) {      else if(b_is_n && (force_num || (maybe_bool && !a_is_n))) {
540          // coerce 'b' to  double          // coerce 'b' to  double
541          SEXP db = PROTECT(nz2Csparse(b, x_double)); nprot++;          SEXP db = PROTECT(nz2Csparse(b, x_double)); nprot++;
542          chb = AS_CHM_SP(db);          chb = AS_CHM_SP(db);
543          R_CheckStack();          R_CheckStack();
544          b_is_n = FALSE;          // b_is_n = FALSE;
545        }
546        else if(do_bool == TRUE) { // Want boolean arithmetic: sufficient if *one* is pattern:
547            if(!a_is_n && !b_is_n) {
548                // coerce 'a' to pattern
549                SEXP da = PROTECT(Csparse2nz(a, /* tri = */
550                                             R_check_class_etc(a, valid_tri) >= 0)); nprot++;
551                cha = AS_CHM_SP(da);
552                R_CheckStack();
553                // a_is_n = TRUE;
554            }
555      }      }
   
556      chTr = cholmod_transpose((tr) ? chb : cha, chb->xtype, &c);      chTr = cholmod_transpose((tr) ? chb : cha, chb->xtype, &c);
557      chc = cholmod_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb,      chc = cholmod_ssmult((tr) ? cha : chTr, (tr) ? chTr : chb,
558                           /*out_stype:*/ 0, /* values : */ do_bool != TRUE,                           /*out_stype:*/ 0, /* values : */ do_bool != TRUE,
# Line 520  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 545  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 589  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 670  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 683  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 691  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 699  Line 742 
742          chx = AS_CHM_SP(dx);          chx = AS_CHM_SP(dx);
743          R_CheckStack();          R_CheckStack();
744      }      }
745        else if(do_bool == TRUE && !x_is_n) { // Want boolean arithmetic; need patter[n]
746            // coerce 'x' to pattern
747            static const char *valid_tri[] = { MATRIX_VALID_tri_Csparse, "" };
748            SEXP dx = PROTECT(Csparse2nz(x, /* tri = */
749                                         R_check_class_etc(x, valid_tri) >= 0)); nprot++;
750            chx = AS_CHM_SP(dx);
751            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 721  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 742  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 774  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 846  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.3069  
changed lines
  Added in v.3147

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business University of Wisconsin - Madison Powered By FusionForge