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 2137, Mon Mar 17 22:21:24 2008 UTC
# Line 475  Line 475 
475      fclose(f);      fclose(f);
476      return R_NilValue;      return R_NilValue;
477  }  }
478    
479    
480    /**
481     * Extract the diagonal entries from *triangular* Csparse matrix  __or__ a
482     * cholmod_sparse factor (LDL = TRUE).
483     *
484     * @param n  dimension of the matrix.
485     * @param x_p  'p' (column pointer) slot contents
486     * @param x_x  'x' (non-zero entries) slot contents
487     * @param perm 'perm' (= permutation vector) slot contents
488     * @param resultKind a (SEXP) string indicating which kind of result is desired.
489     *
490     * @return  a SEXP, either a (double) number or a length n-vector of diagonal entries
491     */
492    SEXP diag_tC_ptr(int n, int *x_p, double *x_x, int *perm, SEXP resultKind)
493    /*                                ^^^^^^ FIXME[Generalize] to int / ... */
494    {
495        const char* res_ch = CHAR(STRING_ELT(resultKind,0));
496        enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log
497        } res_kind = ((!strcmp(res_ch, "trace")) ? trace :
498                      ((!strcmp(res_ch, "sumLog")) ? sum_log :
499                       ((!strcmp(res_ch, "prod")) ? prod :
500                        ((!strcmp(res_ch, "diag")) ? diag :
501                         ((!strcmp(res_ch, "diagBack")) ? diag_backpermuted :
502                          -1)))));
503        int i, n_x, i_from = 0;
504        SEXP ans = PROTECT(allocVector(REALSXP,
505    /*                                 ^^^^  FIXME[Generalize] */
506                                       (res_kind == diag ||
507                                        res_kind == diag_backpermuted) ? n : 1));
508        double *v = REAL(ans);
509    /*  ^^^^^^      ^^^^  FIXME[Generalize] */
510    
511    #define for_DIAG(v_ASSIGN)                                              \
512        for(i = 0; i < n; i++, i_from += n_x) {                             \
513            /* looking at i-th column */                                    \
514            n_x = x_p[i+1] - x_p[i];/* #{entries} in this column */ \
515            v_ASSIGN;                                                       \
516        }
517    
518        /* NOTA BENE: we assume  -- uplo = "L" i.e. lower triangular matrix
519         *            for uplo = "U" (makes sense with a "dtCMatrix" !),
520         *            should use  x_x[i_from + (nx - 1)] instead of x_x[i_from],
521         *            where nx = (x_p[i+1] - x_p[i])
522         */
523    
524        switch(res_kind) {
525        case trace:
526            v[0] = 0.;
527            for_DIAG(v[0] += x_x[i_from]);
528            break;
529    
530        case sum_log:
531            v[0] = 0.;
532            for_DIAG(v[0] += log(x_x[i_from]));
533            break;
534    
535        case prod:
536            v[0] = 1.;
537            for_DIAG(v[0] *= x_x[i_from]);
538            break;
539    
540        case diag:
541            for_DIAG(v[i] = x_x[i_from]);
542            break;
543    
544        case diag_backpermuted:
545            for_DIAG(v[i] = x_x[i_from]);
546            /* now back_permute : */
547    
548            error(_("resultKind = 'diagBack' (back-permuted) is not yet implemented"));
549            break;
550    
551        default: /* -1 from above */
552            error("diag_tC(): invalid 'resultKind'");
553            /* Wall: */ ans = R_NilValue; v = REAL(ans);
554        }
555    
556        UNPROTECT(1);
557        return ans;
558    }
559    
560    /**
561     * Extract the diagonal entries from *triangular* Csparse matrix  __or__ a
562     * cholmod_sparse factor (LDL = TRUE).
563     *
564     * @param pslot  'p' (column pointer)   slot of Csparse matrix/factor
565     * @param xslot  'x' (non-zero entries) slot of Csparse matrix/factor
566     * @param perm_slot  'perm' (= permutation vector) slot of corresponding CHMfactor
567     * @param resultKind a (SEXP) string indicating which kind of result is desired.
568     *
569     * @return  a SEXP, either a (double) number or a length n-vector of diagonal entries
570     */
571    SEXP diag_tC(SEXP pslot, SEXP xslot, SEXP perm_slot, SEXP resultKind)
572    {
573        int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */
574            *x_p  = INTEGER(pslot),
575            *perm = INTEGER(perm_slot);
576        double *x_x = REAL(xslot);
577    /*  ^^^^^^        ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/
578    
579        return diag_tC_ptr(n, x_p, x_x, perm, resultKind);
580    }

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

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