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 2958, Mon Jan 20 16:22:40 2014 UTC revision 2984, Sat Apr 12 21:37:37 2014 UTC
# Line 631  Line 631 
631   *   *
632   * @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
633   */   */
634  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,
635  /*                                ^^^^^^ FIXME[Generalize] to int / ... */  /*                                ^^^^^^ FIXME[Generalize] to int / ... */
636                     SEXP resultKind)
637  {  {
638      const char* res_ch = CHAR(STRING_ELT(resultKind,0));      const char* res_ch = CHAR(STRING_ELT(resultKind,0));
639      enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log      enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log, min, max, range
640      } res_kind = ((!strcmp(res_ch, "trace")) ? trace :      } res_kind = ((!strcmp(res_ch, "trace")) ? trace :
641                    ((!strcmp(res_ch, "sumLog")) ? sum_log :                    ((!strcmp(res_ch, "sumLog")) ? sum_log :
642                     ((!strcmp(res_ch, "prod")) ? prod :                     ((!strcmp(res_ch, "prod")) ? prod :
643                        ((!strcmp(res_ch, "min")) ? min :
644                         ((!strcmp(res_ch, "max")) ? max :
645                          ((!strcmp(res_ch, "range")) ? range :
646                      ((!strcmp(res_ch, "diag")) ? diag :                      ((!strcmp(res_ch, "diag")) ? diag :
647                       ((!strcmp(res_ch, "diagBack")) ? diag_backpermuted :                       ((!strcmp(res_ch, "diagBack")) ? diag_backpermuted :
648                        -1)))));                           -1))))))));
649      int i, n_x, i_from = 0;      int i, n_x, i_from;
650      SEXP ans = PROTECT(allocVector(REALSXP,      SEXP ans = PROTECT(allocVector(REALSXP,
651  /*                                 ^^^^  FIXME[Generalize] */  /*                                 ^^^^  FIXME[Generalize] */
652                                     (res_kind == diag ||                                     (res_kind == diag ||
653                                      res_kind == diag_backpermuted) ? n : 1));                                      res_kind == diag_backpermuted) ? n :
654                                       (res_kind == range ? 2 : 1)));
655      double *v = REAL(ans);      double *v = REAL(ans);
656  /*  ^^^^^^      ^^^^  FIXME[Generalize] */  /*  ^^^^^^      ^^^^  FIXME[Generalize] */
657    
658        i_from = (is_U ? -1 : 0);
659    
660  #define for_DIAG(v_ASSIGN)                                      \  #define for_DIAG(v_ASSIGN)                                      \
661      for(i = 0; i < n; i++, i_from += n_x) {                     \      for(i = 0; i < n; i++) {                                    \
662          /* looking at i-th column */                            \          /* looking at i-th column */                            \
663          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 */ \
664            if( is_U) i_from += n_x;                                \
665          v_ASSIGN;                                               \          v_ASSIGN;                                               \
666            if(!is_U) i_from += n_x;                                \
667      }      }
668    
669      /* NOTA BENE: we assume  -- uplo = "L" i.e. lower triangular matrix      /* NOTA BENE: we assume  -- uplo = "L" i.e. lower triangular matrix
670       *            for uplo = "U" (makes sense with a "dtCMatrix" !),       *            for uplo = "U" (makes sense with a "dtCMatrix" !),
671       *            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],
672       *            where nx = (x_p[i+1] - x_p[i])       *            where n_x = (x_p[i+1] - x_p[i])
673       */       */
674    
675      switch(res_kind) {      switch(res_kind) {
676      case trace:      case trace: // = sum
677          v[0] = 0.;          v[0] = 0.;
678          for_DIAG(v[0] += x_x[i_from]);          for_DIAG(v[0] += x_x[i_from]);
679          break;          break;
# Line 679  Line 688 
688          for_DIAG(v[0] *= x_x[i_from]);          for_DIAG(v[0] *= x_x[i_from]);
689          break;          break;
690    
691        case min:
692            v[0] = R_PosInf;
693            for_DIAG(if(v[0] > x_x[i_from]) v[0] = x_x[i_from]);
694            break;
695    
696        case max:
697            v[0] = R_NegInf;
698            for_DIAG(if(v[0] < x_x[i_from]) v[0] = x_x[i_from]);
699            break;
700    
701        case range:
702            v[0] = R_PosInf;
703            v[1] = R_NegInf;
704            for_DIAG(if(v[0] > x_x[i_from]) v[0] = x_x[i_from];
705                     if(v[1] < x_x[i_from]) v[1] = x_x[i_from]);
706            break;
707    
708      case diag:      case diag:
709          for_DIAG(v[i] = x_x[i_from]);          for_DIAG(v[i] = x_x[i_from]);
710          break;          break;
# Line 708  Line 734 
734   * Extract the diagonal entries from *triangular* Csparse matrix  __or__ a   * Extract the diagonal entries from *triangular* Csparse matrix  __or__ a
735   * cholmod_sparse factor (LDL = TRUE).   * cholmod_sparse factor (LDL = TRUE).
736   *   *
737     * @param obj -- now a cholmod_sparse factor or a dtCMatrix
738   * @param pslot  'p' (column pointer)   slot of Csparse matrix/factor   * @param pslot  'p' (column pointer)   slot of Csparse matrix/factor
739   * @param xslot  'x' (non-zero entries) slot of Csparse matrix/factor   * @param xslot  'x' (non-zero entries) slot of Csparse matrix/factor
740   * @param perm_slot  'perm' (= permutation vector) slot of corresponding CHMfactor;   * @param perm_slot  'perm' (= permutation vector) slot of corresponding CHMfactor;
# Line 716  Line 743 
743   *   *
744   * @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
745   */   */
746  SEXP diag_tC(SEXP pslot, SEXP xslot, SEXP perm_slot, SEXP resultKind)  SEXP diag_tC(SEXP obj, SEXP resultKind)
747  {  {
748    
749        SEXP
750            pslot = GET_SLOT(obj, Matrix_pSym),
751            xslot = GET_SLOT(obj, Matrix_xSym);
752        Rboolean is_U = (R_has_slot(obj, Matrix_uploSym) &&
753                         *CHAR(asChar(GET_SLOT(obj, Matrix_uploSym))) == 'U');
754      int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */      int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */
755          *x_p  = INTEGER(pslot),          *x_p  = INTEGER(pslot), pp = -1, *perm;
         *perm = INTEGER(perm_slot);  
756      double *x_x = REAL(xslot);      double *x_x = REAL(xslot);
757  /*  ^^^^^^        ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/  /*  ^^^^^^        ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/
758    
759      return diag_tC_ptr(n, x_p, x_x, perm, resultKind);      if(R_has_slot(obj, Matrix_permSym))
760            perm = INTEGER(GET_SLOT(obj, Matrix_permSym));
761        else perm = &pp;
762    
763        return diag_tC_ptr(n, x_p, x_x, is_U, perm, resultKind);
764  }  }
765    
766    
767  /**  /**
768   * Create a Csparse matrix object from indices and/or pointers.   * Create a Csparse matrix object from indices and/or pointers.
769   *   *

Legend:
Removed from v.2958  
changed lines
  Added in v.2984

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