SCM

SCM Repository

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

Diff of /pkg/src/Csparse.c

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

revision 2125, Fri Mar 7 07:58:26 2008 UTC revision 2175, Wed Apr 23 11:23:50 2008 UTC
# Line 349  Line 349 
349    
350  SEXP Csparse_drop(SEXP x, SEXP tol)  SEXP Csparse_drop(SEXP x, SEXP tol)
351  {  {
352        const char *cl = class_P(x);
353        /* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */
354        int tr = (cl[1] == 't');
355      CHM_SP chx = AS_CHM_SP(x);      CHM_SP chx = AS_CHM_SP(x);
356      CHM_SP ans = cholmod_copy(chx, chx->stype, chx->xtype, &c);      CHM_SP ans = cholmod_copy(chx, chx->stype, chx->xtype, &c);
357      double dtol = asReal(tol);      double dtol = asReal(tol);
# Line 357  Line 360 
360    
361      if(!cholmod_drop(dtol, ans, &c))      if(!cholmod_drop(dtol, ans, &c))
362          error(_("cholmod_drop() failed"));          error(_("cholmod_drop() failed"));
363      return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "",      return chm_sparse_to_SEXP(ans, 1,
364                                  tr ? ((*uplo_P(x) == 'U') ? 1 : -1) : 0,
365                                  Rkind, tr ? diag_P(x) : "",
366                                GET_SLOT(x, Matrix_DimNamesSym));                                GET_SLOT(x, Matrix_DimNamesSym));
367  }  }
368    
# Line 475  Line 480 
480      fclose(f);      fclose(f);
481      return R_NilValue;      return R_NilValue;
482  }  }
483    
484    
485    /**
486     * Extract the diagonal entries from *triangular* Csparse matrix  __or__ a
487     * cholmod_sparse factor (LDL = TRUE).
488     *
489     * @param n  dimension of the matrix.
490     * @param x_p  'p' (column pointer) slot contents
491     * @param x_x  'x' (non-zero entries) slot contents
492     * @param perm 'perm' (= permutation vector) slot contents; only used for "diagBack"
493     * @param resultKind a (SEXP) string indicating which kind of result is desired.
494     *
495     * @return  a SEXP, either a (double) number or a length n-vector of diagonal entries
496     */
497    SEXP diag_tC_ptr(int n, int *x_p, double *x_x, int *perm, SEXP resultKind)
498    /*                                ^^^^^^ FIXME[Generalize] to int / ... */
499    {
500        const char* res_ch = CHAR(STRING_ELT(resultKind,0));
501        enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log
502        } res_kind = ((!strcmp(res_ch, "trace")) ? trace :
503                      ((!strcmp(res_ch, "sumLog")) ? sum_log :
504                       ((!strcmp(res_ch, "prod")) ? prod :
505                        ((!strcmp(res_ch, "diag")) ? diag :
506                         ((!strcmp(res_ch, "diagBack")) ? diag_backpermuted :
507                          -1)))));
508        int i, n_x, i_from = 0;
509        SEXP ans = PROTECT(allocVector(REALSXP,
510    /*                                 ^^^^  FIXME[Generalize] */
511                                       (res_kind == diag ||
512                                        res_kind == diag_backpermuted) ? n : 1));
513        double *v = REAL(ans);
514    /*  ^^^^^^      ^^^^  FIXME[Generalize] */
515    
516    #define for_DIAG(v_ASSIGN)                                              \
517        for(i = 0; i < n; i++, i_from += n_x) {                             \
518            /* looking at i-th column */                                    \
519            n_x = x_p[i+1] - x_p[i];/* #{entries} in this column */ \
520            v_ASSIGN;                                                       \
521        }
522    
523        /* NOTA BENE: we assume  -- uplo = "L" i.e. lower triangular matrix
524         *            for uplo = "U" (makes sense with a "dtCMatrix" !),
525         *            should use  x_x[i_from + (nx - 1)] instead of x_x[i_from],
526         *            where nx = (x_p[i+1] - x_p[i])
527         */
528    
529        switch(res_kind) {
530        case trace:
531            v[0] = 0.;
532            for_DIAG(v[0] += x_x[i_from]);
533            break;
534    
535        case sum_log:
536            v[0] = 0.;
537            for_DIAG(v[0] += log(x_x[i_from]));
538            break;
539    
540        case prod:
541            v[0] = 1.;
542            for_DIAG(v[0] *= x_x[i_from]);
543            break;
544    
545        case diag:
546            for_DIAG(v[i] = x_x[i_from]);
547            break;
548    
549        case diag_backpermuted:
550            for_DIAG(v[i] = x_x[i_from]);
551    
552            warning(_("resultKind = 'diagBack' (back-permuted) is experimental"));
553            /* now back_permute : */
554            for(i = 0; i < n; i++) {
555                double tmp = v[i]; v[i] = v[perm[i]]; v[perm[i]] = tmp;
556                /*^^^^ FIXME[Generalize] */
557            }
558            break;
559    
560        default: /* -1 from above */
561            error("diag_tC(): invalid 'resultKind'");
562            /* Wall: */ ans = R_NilValue; v = REAL(ans);
563        }
564    
565        UNPROTECT(1);
566        return ans;
567    }
568    
569    /**
570     * Extract the diagonal entries from *triangular* Csparse matrix  __or__ a
571     * cholmod_sparse factor (LDL = TRUE).
572     *
573     * @param pslot  'p' (column pointer)   slot of Csparse matrix/factor
574     * @param xslot  'x' (non-zero entries) slot of Csparse matrix/factor
575     * @param perm_slot  'perm' (= permutation vector) slot of corresponding CHMfactor;
576     *                   only used for "diagBack"
577     * @param resultKind a (SEXP) string indicating which kind of result is desired.
578     *
579     * @return  a SEXP, either a (double) number or a length n-vector of diagonal entries
580     */
581    SEXP diag_tC(SEXP pslot, SEXP xslot, SEXP perm_slot, SEXP resultKind)
582    {
583        int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */
584            *x_p  = INTEGER(pslot),
585            *perm = INTEGER(perm_slot);
586        double *x_x = REAL(xslot);
587    /*  ^^^^^^        ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/
588    
589        return diag_tC_ptr(n, x_p, x_x, perm, resultKind);
590    }

Legend:
Removed from v.2125  
changed lines
  Added in v.2175

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