SCM

SCM Repository

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

Diff of /pkg/src/bCrosstab.c

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

revision 356, Sat Nov 27 01:45:39 2004 UTC revision 359, Sat Nov 27 17:58:23 2004 UTC
# Line 1  Line 1 
 /* TODO: 1) combine diag_update and offdiag_update - DONE  
          2) Move Parent_inverse from tscMatrix to here  
          3) Modify the creation of the L diagonals to include a call  
          to ldl_numeric to get Li.  
 */  
   
1  #include "bCrosstab.h"  #include "bCrosstab.h"
2    
3  /**  /**
# Line 332  Line 326 
326          int dind = Lind(j, j), i, k;          int dind = Lind(j, j), i, k;
327          SEXP ctd = VECTOR_ELT(ctab, dind); /* diagonal in crosstab */          SEXP ctd = VECTOR_ELT(ctab, dind); /* diagonal in crosstab */
328          SEXP Dimslot = GET_SLOT(ctd, Matrix_DimSym),          SEXP Dimslot = GET_SLOT(ctd, Matrix_DimSym),
329              Ljj, cpp = GET_SLOT(ctd, Matrix_pSym),              Linvj, Ljj, cpp = GET_SLOT(ctd, Matrix_pSym),
330              cip = GET_SLOT(ctd, Matrix_iSym);              cip = GET_SLOT(ctd, Matrix_iSym);
331          int *Lp, *Perm, *cp = INTEGER(cpp),          int *Lp, *Linvp, *Perm, *cp = INTEGER(cpp),
332              *ci = INTEGER(cip), *parent,              *ci = INTEGER(cip), *parent,
333              ncj = length(cpp) - 1,              ncj = length(cpp) - 1,
334              nnz = length(cip);              nnz = length(cip);
335            double *dtmp;
336    
337          SET_VECTOR_ELT(Parent, j, allocVector(INTSXP, ncj));          SET_VECTOR_ELT(Parent, j, allocVector(INTSXP, ncj));
338          parent = INTEGER(VECTOR_ELT(Parent, j));          parent = INTEGER(VECTOR_ELT(Parent, j));
# Line 345  Line 340 
340          Perm = INTEGER(VECTOR_ELT(perm, j));          Perm = INTEGER(VECTOR_ELT(perm, j));
341          SET_VECTOR_ELT(L, dind, NEW_OBJECT(MAKE_CLASS("cscMatrix")));          SET_VECTOR_ELT(L, dind, NEW_OBJECT(MAKE_CLASS("cscMatrix")));
342          Ljj = VECTOR_ELT(L, dind);          Ljj = VECTOR_ELT(L, dind);
343            SET_VECTOR_ELT(Linv, j, NEW_OBJECT(MAKE_CLASS("cscMatrix")));
344            Linvj = VECTOR_ELT(Linv, j);
345          SET_SLOT(Ljj, Matrix_DimSym, duplicate(Dimslot));          SET_SLOT(Ljj, Matrix_DimSym, duplicate(Dimslot));
346            SET_SLOT(Linvj, Matrix_DimSym, duplicate(Dimslot));
347          SET_SLOT(Ljj, Matrix_factorization, allocVector(VECSXP, 0));          SET_SLOT(Ljj, Matrix_factorization, allocVector(VECSXP, 0));
348            SET_SLOT(Linvj, Matrix_factorization, allocVector(VECSXP, 0));
349          SET_SLOT(Ljj, Matrix_pSym, allocVector(INTSXP, ncj + 1));          SET_SLOT(Ljj, Matrix_pSym, allocVector(INTSXP, ncj + 1));
350            SET_SLOT(Linvj, Matrix_pSym, allocVector(INTSXP, ncj + 1));
351          Lp = INTEGER(GET_SLOT(Ljj, Matrix_pSym));          Lp = INTEGER(GET_SLOT(Ljj, Matrix_pSym));
352            Linvp = INTEGER(GET_SLOT(Linvj, Matrix_pSym));
353          if (nnz > ncj) {        /* calculate fill-reducing permutation */          if (nnz > ncj) {        /* calculate fill-reducing permutation */
354              int *iPerm = Calloc(ncj, int),              int *Li, *Lnz = Calloc(ncj, int), *tmp = Calloc(ncj, int), info;
                 *Lnz = Calloc(ncj, int);  
             double *Lx;  
355    
356              ssc_metis_order(ncj, cp, ci, Perm, iPerm);              ssc_metis_order(ncj, cp, ci, Perm, tmp);
357                                  /* apply to the crosstabulation */                                  /* apply to the crosstabulation */
358              bCrosstab_permute(ctab, nf, j, Perm, iPerm);              bCrosstab_permute(ctab, nf, j, Perm, tmp);
359                                  /* apply to the factor */                                  /* apply to the factor */
360              factor_levels_permute(VECTOR_ELT(fcp, j),              factor_levels_permute(VECTOR_ELT(fcp, j),
361                                    VECTOR_ELT(flist, j), Perm, iPerm);                                    VECTOR_ELT(flist, j), Perm, tmp);
362                                  /* symbolic analysis to get Parent */                                  /* symbolic analysis to get Parent */
363              ldl_symbolic(ncj, cp, ci, Lp, /* iPerm used as scratch */              ldl_symbolic(ncj, cp, ci, Lp, parent, Lnz, tmp,
                          parent, Lnz, iPerm,  
364                           (int *) NULL, (int *) NULL);                           (int *) NULL, (int *) NULL);
365                                    /* decompose the identity to get the row pointers */
366                dtmp = REAL(GET_SLOT(ctd, Matrix_xSym));
367                for (i = 0; i < nnz; i++) dtmp[i] = 0.; /* initialize */
368                                    /* diagonal el is last element in the column */
369                for (i = 0; i < ncj; i++) dtmp[cp[i+1] - 1] = 1.;
370              nnz = Lp[ncj];              nnz = Lp[ncj];
             Free(iPerm); Free(Lnz);  
371              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, nnz));              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, nnz));
372              SET_SLOT(Ljj, Matrix_xSym, allocVector(REALSXP, nnz));              SET_SLOT(Ljj, Matrix_xSym, allocVector(REALSXP, nnz));
373              Lnz = INTEGER(GET_SLOT(Ljj, Matrix_iSym));              info = ldl_numeric(ncj, cp, ci, dtmp, Lp, parent, Lnz,
374              Lx = REAL(GET_SLOT(Ljj, Matrix_xSym));                                 INTEGER(GET_SLOT(Ljj, Matrix_iSym)),
375              for (i = 0; i < nnz; i++) {                                 REAL(GET_SLOT(Ljj, Matrix_xSym)),
376                  Lnz[i] = 0; Lx[i] = 0.;                                 (double *) R_alloc(ncj, sizeof(double)), /* D */
377              }                                 (double *) R_alloc(ncj, sizeof(double)), /* Y */
378                                   (int *) R_alloc(ncj, sizeof(int)),       /* Pattern */
379                                   (int *) R_alloc(ncj, sizeof(int)),       /* Flag */
380                                   (int *) NULL, (int *) NULL);
381                if (info < ncj) error("identity matrix is not positive definite?");
382                Free(Lnz); Free(tmp);
383          } else {          } else {
384              for (i = 0; i <= ncj; i++) {              for (i = 0; i <= ncj; i++) {
385                  Lp[i] = 0;                  Lp[i] = 0;
# Line 382  Line 389 
389              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, 0));              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, 0));
390              SET_SLOT(Ljj, Matrix_xSym, allocVector(REALSXP, 0));              SET_SLOT(Ljj, Matrix_xSym, allocVector(REALSXP, 0));
391          }          }
392          SET_VECTOR_ELT(Linv, j,                                  /* Derive the diagonal block of Linv */
393                         Parent_inverse(VECTOR_ELT(Parent, j),          nnz = parent_inv_ap(ncj, 0, parent, Linvp);
394                                        ScalarLogical(1)));          SET_SLOT(Linvj, Matrix_iSym, allocVector(INTSXP, nnz));
395            parent_inv_ai(ncj, 0, parent, INTEGER(GET_SLOT(Linvj, Matrix_iSym)));
396            SET_SLOT(Linvj, Matrix_xSym, allocVector(REALSXP, nnz));
397            dtmp = REAL(GET_SLOT(Linvj, Matrix_xSym));
398            for (i = 0; i < nnz; i++) dtmp[i] = 1.;
399  /* FIXME: Update any blocks below the diagonal in this column if necessary*/  /* FIXME: Update any blocks below the diagonal in this column if necessary*/
400          for (k = j+1; k < nf; k++) { /* update remaining columns, if any */          for (k = j+1; k < nf; k++) { /* update remaining columns, if any */
401              for (i = k; i < nf; i++) block_update(ctab, j, k, i);              for (i = k; i < nf; i++) block_update(ctab, j, k, i);

Legend:
Removed from v.356  
changed lines
  Added in v.359

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