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 3072, Fri Mar 27 15:10:48 2015 UTC revision 3076, Mon Mar 30 10:23:42 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 409  Line 412 
412                                Rkind, tr ? diag_P(x) : "", dn);                                Rkind, tr ? diag_P(x) : "", dn);
413  }  }
414    
415  /* NOTA BENE:  cholmod_ssmult(A,B, ...) ->  ./CHOLMOD/MatrixOps/cholmod_ssmult.c  /** @brief  A %*% B  - for matrices of class CsparseMatrix (R package "Matrix")
416     *
417     * @param a
418     * @param b
419     * @param bool_arith
420     *
421     * @return
422     *
423     * NOTA BENE:  cholmod_ssmult(A,B, ...) ->  ./CHOLMOD/MatrixOps/cholmod_ssmult.c
424   * ---------  computes a patter*n* matrix __always_ when   * ---------  computes a patter*n* matrix __always_ when
425   * *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):
426     ---------------------------------------------------------------------------     ---------------------------------------------------------------------------
# Line 418  Line 429 
429     ---------------------------------------------------------------------------     ---------------------------------------------------------------------------
430   * ==> 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 !!!
431   */   */
   
432  SEXP Csparse_Csparse_prod(SEXP a, SEXP b, SEXP bool_arith)  SEXP Csparse_Csparse_prod(SEXP a, SEXP b, SEXP bool_arith)
433  {  {
434      CHM_SP      CHM_SP
# Line 484  Line 494 
494      return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn);      return chm_sparse_to_SEXP(chc, 1, uploT, /*Rkind*/0, diag, dn);
495  }  }
496    
497  /* trans = FALSE:  crossprod(a,b)  /** @brief [t]crossprod (<Csparse>, <Csparse>)
498   * trans = TRUE : tcrossprod(a,b) */   *
499     * @param a a "CsparseMatrix" object
500     * @param b a "CsparseMatrix" object
501     * @param trans trans = FALSE:  crossprod(a,b)
502     *              trans = TRUE : tcrossprod(a,b)
503     * @param bool_arith logical (TRUE / NA / FALSE): Should boolean arithmetic be used.
504     *
505     * @return a CsparseMatrix, the (t)cross product of a and b.
506     */
507  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)
508  {  {
509      int tr = asLogical(trans), nprot = 1,      int tr = asLogical(trans), nprot = 1,
# Line 561  Line 579 
579  /**  /**
580   * All (dense * sparse)  Matrix products and cross products   * All (dense * sparse)  Matrix products and cross products
581   *   *
582   *   f( f(<Csparse>)  %*%  f(<dense>) )   for  f in {t(), identity()}   *   f( f(<Csparse>)  %*%  f(<dense>) )   where  f ()  is either t () [tranpose] or the identity.
583   *   *
584   * @param a CsparseMatrix  (n x m)   * @param a CsparseMatrix  (n x m)
585   * @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 686  Line 704 
704  }  }
705    
706    
707  /* Computes   x'x  or  x x' -- *also* for Tsparse (triplet = TRUE)  /** @brief Computes   x'x  or  x x' -- *also* for Tsparse (triplet = TRUE)
708     see Csparse_Csparse_crossprod above for  x'y and x y' */      see Csparse_Csparse_crossprod above for  x'y and x y'
709    */
710  SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet, SEXP bool_arith)  SEXP Csparse_crossprod(SEXP x, SEXP trans, SEXP triplet, SEXP bool_arith)
711  {  {
712      int tripl = asLogical(triplet),      int tripl = asLogical(triplet),
# Line 699  Line 718 
718      SEXP xx = PROTECT(Tsparse_diagU2N(x));      SEXP xx = PROTECT(Tsparse_diagU2N(x));
719      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;
720  #endif  #endif
721      CHM_SP chcp, chxt,      CHM_SP chcp, chxt, chxc,
722          chx = (tripl ?          chx = (tripl ?
723                 cholmod_triplet_to_sparse(cht, cht->nnz, &c) :                 cholmod_triplet_to_sparse(cht, cht->nnz, &c) :
724                 AS_CHM_SP(x));                 AS_CHM_SP(x));
# Line 707  Line 726 
726      R_CheckStack();      R_CheckStack();
727      Rboolean      Rboolean
728          x_is_n = (chx->xtype == CHOLMOD_PATTERN),          x_is_n = (chx->xtype == CHOLMOD_PATTERN),
729            x_is_sym = chx->stype != 0,
730          force_num = (do_bool == FALSE);          force_num = (do_bool == FALSE);
731    
732      if(x_is_n && force_num) {      if(x_is_n && force_num) {
# Line 726  Line 746 
746    
747      if (!tr) chxt = cholmod_transpose(chx, chx->xtype, &c);      if (!tr) chxt = cholmod_transpose(chx, chx->xtype, &c);
748    
749      chcp = cholmod_aat((!tr) ? chxt : chx, (int *) NULL, 0,      if (x_is_sym) // cholmod_aat() does not like symmetric
750                         /* mode: */ chx->xtype, &c);          chxc = cholmod_copy(tr ? chx : chxt, /* stype: */ 0,
751                                chx->xtype, &c);
752        // CHOLMOD/Core/cholmod_aat.c :
753        chcp = cholmod_aat(x_is_sym ? chxc : (tr ? chx : chxt),
754                           (int *) NULL, 0, /* mode: */ chx->xtype, &c);
755      if(!chcp) {      if(!chcp) {
756          UNPROTECT(1);          UNPROTECT(1);
757          error(_("Csparse_crossprod(): error return from cholmod_aat()"));          error(_("Csparse_crossprod(): error return from cholmod_aat()"));
# Line 745  Line 769 
769      return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn);      return chm_sparse_to_SEXP(chcp, 1, 0, 0, "", dn);
770  }  }
771    
772  /* 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,
773  *  at least all "explicit" zeros */   *  at least all "explicit" zeros. */
774  SEXP Csparse_drop(SEXP x, SEXP tol)  SEXP Csparse_drop(SEXP x, SEXP tol)
775  {  {
776      const char *cl = class_P(x);      const char *cl = class_P(x);
# Line 766  Line 790 
790                                GET_SLOT(x, Matrix_DimNamesSym));                                GET_SLOT(x, Matrix_DimNamesSym));
791  }  }
792    
793    /** @brief Horizontal Concatenation -  cbind( <Csparse>,  <Csparse>)
794     */
795  SEXP Csparse_horzcat(SEXP x, SEXP y)  SEXP Csparse_horzcat(SEXP x, SEXP y)
796  {  {
797  #define CSPARSE_CAT(_KIND_)                                             \  #define CSPARSE_CAT(_KIND_)                                             \
# Line 798  Line 824 
824                                1, 0, Rkind, "", R_NilValue);                                1, 0, Rkind, "", R_NilValue);
825  }  }
826    
827    /** @brief Vertical Concatenation -  rbind( <Csparse>,  <Csparse>)
828     */
829  SEXP Csparse_vertcat(SEXP x, SEXP y)  SEXP Csparse_vertcat(SEXP x, SEXP y)
830  {  {
831      CSPARSE_CAT("vertcat");      CSPARSE_CAT("vertcat");
# Line 870  Line 898 
898  }  }
899    
900  /**  /**
901   * "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
902   * Working via CHOLMOD_submatrix, see ./CHOLMOD/MatrixOps/cholmod_submatrix.c   * Working via CHOLMOD_submatrix, see ./CHOLMOD/MatrixOps/cholmod_submatrix.c
903   * @param x CsparseMatrix   * @param x CsparseMatrix
904   * @param i row     indices (0-origin), or NULL (R's)   * @param i row     indices (0-origin), or NULL (R, not C)
905   * @param j columns indices (0-origin), or NULL   * @param j columns indices (0-origin), or NULL
906   *   *
907   * @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.3076

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