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 2901, Tue Sep 10 10:45:25 2013 UTC revision 3018, Sat Oct 11 17:52:10 2014 UTC
# Line 186  Line 186 
186      return ans;      return ans;
187  }  }
188    
189  SEXP Csparse_to_matrix(SEXP x)  SEXP Csparse_to_matrix(SEXP x, SEXP chk)
190  {  {
191      return chm_dense_to_matrix(cholmod_sparse_to_dense(AS_CHM_SP__(x), &c),      return chm_dense_to_matrix(cholmod_sparse_to_dense(AS_CHM_SP2(x, asLogical(chk)), &c),
192                                 1 /*do_free*/, GET_SLOT(x, Matrix_DimNamesSym));                                 1 /*do_free*/, GET_SLOT(x, Matrix_DimNamesSym));
193  }  }
194  SEXP Csparse_to_vector(SEXP x)  SEXP Csparse_to_vector(SEXP x)
# Line 345  Line 345 
345  SEXP Csparse_dense_prod(SEXP a, SEXP b)  SEXP Csparse_dense_prod(SEXP a, SEXP b)
346  {  {
347      CHM_SP cha = AS_CHM_SP(a);      CHM_SP cha = AS_CHM_SP(a);
348      SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b));      SEXP b_M = PROTECT(mMatrix_as_dgeMatrix2(b, // transpose_if_vector =
349                                                 cha->ncol == 1));
350      CHM_DN chb = AS_CHM_DN(b_M);      CHM_DN chb = AS_CHM_DN(b_M);
351      CHM_DN chc = cholmod_allocate_dense(cha->nrow, chb->ncol, cha->nrow,      CHM_DN chc = cholmod_allocate_dense(cha->nrow, chb->ncol, cha->nrow,
352                                          chb->xtype, &c);                                          chb->xtype, &c);
# Line 376  Line 377 
377  SEXP Csparse_dense_crossprod(SEXP a, SEXP b)  SEXP Csparse_dense_crossprod(SEXP a, SEXP b)
378  {  {
379      CHM_SP cha = AS_CHM_SP(a);      CHM_SP cha = AS_CHM_SP(a);
380      SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b));      SEXP b_M = PROTECT(mMatrix_as_dgeMatrix2(b, // transpose_if_vector =
381                                                 cha->nrow == 1));
382      CHM_DN chb = AS_CHM_DN(b_M);      CHM_DN chb = AS_CHM_DN(b_M);
383      CHM_DN chc = cholmod_allocate_dense(cha->ncol, chb->ncol, cha->ncol,      CHM_DN chc = cholmod_allocate_dense(cha->ncol, chb->ncol, cha->ncol,
384                                          chb->xtype, &c);                                          chb->xtype, &c);
# Line 569  Line 571 
571      if (csize >= 0 && !isInteger(j))      if (csize >= 0 && !isInteger(j))
572          error(_("Index j must be NULL or integer"));          error(_("Index j must be NULL or integer"));
573    
574    #define CHM_SUB(_M_, _i_, _j_)                                  \
575        cholmod_submatrix(_M_,                                      \
576                          (rsize < 0) ? NULL : INTEGER(_i_), rsize, \
577                          (csize < 0) ? NULL : INTEGER(_j_), csize, \
578                          TRUE, TRUE, &c)
579        CHM_SP ans;
580      if (!chx->stype) {/* non-symmetric Matrix */      if (!chx->stype) {/* non-symmetric Matrix */
581          return chm_sparse_to_SEXP(cholmod_submatrix(chx,          ans = CHM_SUB(chx, i, j);
                                                     (rsize < 0) ? NULL : INTEGER(i), rsize,  
                                                     (csize < 0) ? NULL : INTEGER(j), csize,  
                                                     TRUE, TRUE, &c),  
                                   1, 0, Rkind, "",  
                                   /* FIXME: drops dimnames */ R_NilValue);  
582      }      }
583        else {
584                                  /* for now, cholmod_submatrix() only accepts "generalMatrix" */                                  /* for now, cholmod_submatrix() only accepts "generalMatrix" */
585      CHM_SP tmp = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);      CHM_SP tmp = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);
586      CHM_SP ans = cholmod_submatrix(tmp,          ans = CHM_SUB(tmp, i, j);
                                    (rsize < 0) ? NULL : INTEGER(i), rsize,  
                                    (csize < 0) ? NULL : INTEGER(j), csize,  
                                    TRUE, TRUE, &c);  
587      cholmod_free_sparse(&tmp, &c);      cholmod_free_sparse(&tmp, &c);
588      return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue);      }
589    
590        // "FIXME": currently dropping dimnames, and adding them afterwards in R :
591        /* // dimnames: */
592        /* SEXP x_dns = GET_SLOT(x, Matrix_DimNamesSym), */
593        /*  dn = PROTECT(allocVector(VECSXP, 2)); */
594        return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", /* dimnames: */ R_NilValue);
595  }  }
596    
597  #define _d_Csp_  #define _d_Csp_
# Line 631  Line 638 
638   *   *
639   * @return  a SEXP, either a (double) number or a length n-vector of diagonal entries   * @return  a SEXP, either a (double) number or a length n-vector of diagonal entries
640   */   */
641  SEXP diag_tC_ptr(int n, int *x_p, double *x_x, int *perm, SEXP resultKind)  SEXP diag_tC_ptr(int n, int *x_p, double *x_x, Rboolean is_U, int *perm,
642  /*                                ^^^^^^ FIXME[Generalize] to int / ... */  /*                                ^^^^^^ FIXME[Generalize] to int / ... */
643                     SEXP resultKind)
644  {  {
645      const char* res_ch = CHAR(STRING_ELT(resultKind,0));      const char* res_ch = CHAR(STRING_ELT(resultKind,0));
646      enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log      enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log, min, max, range
647      } res_kind = ((!strcmp(res_ch, "trace")) ? trace :      } res_kind = ((!strcmp(res_ch, "trace")) ? trace :
648                    ((!strcmp(res_ch, "sumLog")) ? sum_log :                    ((!strcmp(res_ch, "sumLog")) ? sum_log :
649                     ((!strcmp(res_ch, "prod")) ? prod :                     ((!strcmp(res_ch, "prod")) ? prod :
650                        ((!strcmp(res_ch, "min")) ? min :
651                         ((!strcmp(res_ch, "max")) ? max :
652                          ((!strcmp(res_ch, "range")) ? range :
653                      ((!strcmp(res_ch, "diag")) ? diag :                      ((!strcmp(res_ch, "diag")) ? diag :
654                       ((!strcmp(res_ch, "diagBack")) ? diag_backpermuted :                       ((!strcmp(res_ch, "diagBack")) ? diag_backpermuted :
655                        -1)))));                           -1))))))));
656      int i, n_x, i_from = 0;      int i, n_x, i_from;
657      SEXP ans = PROTECT(allocVector(REALSXP,      SEXP ans = PROTECT(allocVector(REALSXP,
658  /*                                 ^^^^  FIXME[Generalize] */  /*                                 ^^^^  FIXME[Generalize] */
659                                     (res_kind == diag ||                                     (res_kind == diag ||
660                                      res_kind == diag_backpermuted) ? n : 1));                                      res_kind == diag_backpermuted) ? n :
661                                       (res_kind == range ? 2 : 1)));
662      double *v = REAL(ans);      double *v = REAL(ans);
663  /*  ^^^^^^      ^^^^  FIXME[Generalize] */  /*  ^^^^^^      ^^^^  FIXME[Generalize] */
664    
665        i_from = (is_U ? -1 : 0);
666    
667  #define for_DIAG(v_ASSIGN)                                      \  #define for_DIAG(v_ASSIGN)                                      \
668      for(i = 0; i < n; i++, i_from += n_x) {                     \      for(i = 0; i < n; i++) {                                    \
669          /* looking at i-th column */                            \          /* looking at i-th column */                            \
670          n_x = x_p[i+1] - x_p[i];/* #{entries} in this column */ \          n_x = x_p[i+1] - x_p[i];/* #{entries} in this column */ \
671            if( is_U) i_from += n_x;                                \
672          v_ASSIGN;                                               \          v_ASSIGN;                                               \
673            if(!is_U) i_from += n_x;                                \
674      }      }
675    
676      /* NOTA BENE: we assume  -- uplo = "L" i.e. lower triangular matrix      /* NOTA BENE: we assume  -- uplo = "L" i.e. lower triangular matrix
677       *            for uplo = "U" (makes sense with a "dtCMatrix" !),       *            for uplo = "U" (makes sense with a "dtCMatrix" !),
678       *            should use  x_x[i_from + (nx - 1)] instead of x_x[i_from],       *            should use  x_x[i_from + (n_x - 1)] instead of x_x[i_from],
679       *            where nx = (x_p[i+1] - x_p[i])       *            where n_x = (x_p[i+1] - x_p[i])
680       */       */
681    
682      switch(res_kind) {      switch(res_kind) {
683      case trace:      case trace: // = sum
684          v[0] = 0.;          v[0] = 0.;
685          for_DIAG(v[0] += x_x[i_from]);          for_DIAG(v[0] += x_x[i_from]);
686          break;          break;
# Line 679  Line 695 
695          for_DIAG(v[0] *= x_x[i_from]);          for_DIAG(v[0] *= x_x[i_from]);
696          break;          break;
697    
698        case min:
699            v[0] = R_PosInf;
700            for_DIAG(if(v[0] > x_x[i_from]) v[0] = x_x[i_from]);
701            break;
702    
703        case max:
704            v[0] = R_NegInf;
705            for_DIAG(if(v[0] < x_x[i_from]) v[0] = x_x[i_from]);
706            break;
707    
708        case range:
709            v[0] = R_PosInf;
710            v[1] = R_NegInf;
711            for_DIAG(if(v[0] > x_x[i_from]) v[0] = x_x[i_from];
712                     if(v[1] < x_x[i_from]) v[1] = x_x[i_from]);
713            break;
714    
715      case diag:      case diag:
716          for_DIAG(v[i] = x_x[i_from]);          for_DIAG(v[i] = x_x[i_from]);
717          break;          break;
# Line 708  Line 741 
741   * Extract the diagonal entries from *triangular* Csparse matrix  __or__ a   * Extract the diagonal entries from *triangular* Csparse matrix  __or__ a
742   * cholmod_sparse factor (LDL = TRUE).   * cholmod_sparse factor (LDL = TRUE).
743   *   *
744     * @param obj -- now a cholmod_sparse factor or a dtCMatrix
745   * @param pslot  'p' (column pointer)   slot of Csparse matrix/factor   * @param pslot  'p' (column pointer)   slot of Csparse matrix/factor
746   * @param xslot  'x' (non-zero entries) slot of Csparse matrix/factor   * @param xslot  'x' (non-zero entries) slot of Csparse matrix/factor
747   * @param perm_slot  'perm' (= permutation vector) slot of corresponding CHMfactor;   * @param perm_slot  'perm' (= permutation vector) slot of corresponding CHMfactor;
# Line 716  Line 750 
750   *   *
751   * @return  a SEXP, either a (double) number or a length n-vector of diagonal entries   * @return  a SEXP, either a (double) number or a length n-vector of diagonal entries
752   */   */
753  SEXP diag_tC(SEXP pslot, SEXP xslot, SEXP perm_slot, SEXP resultKind)  SEXP diag_tC(SEXP obj, SEXP resultKind)
754  {  {
755    
756        SEXP
757            pslot = GET_SLOT(obj, Matrix_pSym),
758            xslot = GET_SLOT(obj, Matrix_xSym);
759        Rboolean is_U = (R_has_slot(obj, Matrix_uploSym) &&
760                         *CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))) == 'U');
761      int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */      int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */
762          *x_p  = INTEGER(pslot),          *x_p  = INTEGER(pslot), pp = -1, *perm;
         *perm = INTEGER(perm_slot);  
763      double *x_x = REAL(xslot);      double *x_x = REAL(xslot);
764  /*  ^^^^^^        ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/  /*  ^^^^^^        ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/
765    
766      return diag_tC_ptr(n, x_p, x_x, perm, resultKind);      if(R_has_slot(obj, Matrix_permSym))
767            perm = INTEGER(GET_SLOT(obj, Matrix_permSym));
768        else perm = &pp;
769    
770        return diag_tC_ptr(n, x_p, x_x, is_U, perm, resultKind);
771  }  }
772    
773    
774  /**  /**
775   * Create a Csparse matrix object from indices and/or pointers.   * Create a Csparse matrix object from indices and/or pointers.
776   *   *

Legend:
Removed from v.2901  
changed lines
  Added in v.3018

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