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 444, Wed Jan 19 19:05:15 2005 UTC revision 447, Fri Jan 21 12:17:56 2005 UTC
# Line 293  Line 293 
293   * @param i row index of block to be updated (j < k <= i)   * @param i row index of block to be updated (j < k <= i)
294   */   */
295  static void  static void
296  block_update(SEXP ctab, int j, int k, int i)  block_update(SEXP L, SEXP ZZpO, int j, int k, int i)
297  {  {
298      SEXP tb = VECTOR_ELT(ctab, Lind(i, k)),      SEXP tb = (i == k) ? VECTOR_ELT(ZZpO, i) : VECTOR_ELT(L, Lind(i, k)),
299          ib = VECTOR_ELT(ctab, Lind(i, j)),          ib = VECTOR_ELT(L, Lind(i, j)),
300          kb = VECTOR_ELT(ctab, Lind(k, j));          kb = VECTOR_ELT(L, Lind(k, j));
301      SEXP tpp = GET_SLOT(tb, Matrix_pSym),      SEXP tpp = GET_SLOT(tb, Matrix_pSym),
302          kpp = GET_SLOT(kb, Matrix_pSym);          kpp = GET_SLOT(kb, Matrix_pSym);
303      int *ti = INTEGER(GET_SLOT(tb, Matrix_iSym)),      int *ti = INTEGER(GET_SLOT(tb, Matrix_iSym)),
# Line 384  Line 384 
384   * @param pperm inverse of the permutation   * @param pperm inverse of the permutation
385   */   */
386  static void  static void
387  bCrosstab_permute(SEXP ctab, int nf, int jj, const int nlev[],  bCrosstab_permute(SEXP ctab, int nf, int jj, const int nlev[], const int iperm[])
                   const int perm[], const int iperm[])  
388  {  {
389      int j;      int j;
390      for (j = 0; j < nf; j++) {      for (j = 0; j < nf; j++) {
# Line 409  Line 408 
408      }      }
409  }  }
410    
411    static void
412    symmetric_permute(SEXP A, int nlev, const int iperm[])
413    {
414        SEXP AiP = GET_SLOT(A, Matrix_iSym);
415        int *Ap = INTEGER(GET_SLOT(A, Matrix_pSym)),
416            nnz = length(AiP);
417        double *Ax = REAL(GET_SLOT(A, Matrix_xSym));
418        int *mj = expand_column_pointers(nlev, Ap, Calloc(nnz, int));
419        int *mi = Memcpy(Calloc(nnz, int), INTEGER(AiP), nnz);
420        double *mx = Memcpy(Calloc(nnz, double), Ax, nnz);
421    
422        ind_permute(mi, nnz, iperm);
423        ind_permute(mj, nnz, iperm);
424        make_upper_triangular(mi, mj, nnz);
425        triplet_to_col(nlev, nlev, nnz, mi, mj, mx, Ap, INTEGER(AiP), Ax);
426        Free(mi); Free(mj); Free(mx);
427    }
428    
429  /**  /**
430   * Apply a permutation vector to the levels of a factor.   * Apply a permutation vector to the levels of a factor.
431   *   *
# Line 523  Line 540 
540              int *iPerm = Calloc(ncj, int);              int *iPerm = Calloc(ncj, int);
541    
542              ssc_metis_order(ncj, cp, ci, Perm, iPerm);              ssc_metis_order(ncj, cp, ci, Perm, iPerm);
543                                  /* apply to the crosstabulation and ZZpO */                                  /* apply to the crosstabulation, L, and ZZpO */
544              bCrosstab_permute(ZtZ, nf, j, nlev, Perm, iPerm);              bCrosstab_permute(ZtZ, nf, j, nlev, iPerm);
545              bCrosstab_permute(ZZpO, nf, j, nlev, Perm, iPerm);              bCrosstab_permute(L, nf, j, nlev, iPerm);
546                symmetric_permute(VECTOR_ELT(ZZpO, j), nlev[j], iPerm);
547                                  /* apply to the factor */                                  /* apply to the factor */
548              factor_levels_permute(fac, fcp, Perm, iPerm);              factor_levels_permute(fac, fcp, Perm, iPerm);
549                                  /* symbolic analysis to get Parent */                                  /* symbolic analysis to get Parent */
# Line 554  Line 572 
572                                           VECTOR_ELT(L, Lind(k,j)));                                           VECTOR_ELT(L, Lind(k,j)));
573          }          }
574          for (k = j+1; k < nf; k++) { /* Update remaining columns */          for (k = j+1; k < nf; k++) { /* Update remaining columns */
575              for (i = k; i < nf; i++) block_update(L, j, k, i);              for (i = k; i < nf; i++) block_update(L, ZZpO, j, k, i);
576          }          }
577      }      }
578      /* Convert blockwise Parent arrays to extended Parent arrays */      /* Convert blockwise Parent arrays to extended Parent arrays */

Legend:
Removed from v.444  
changed lines
  Added in v.447

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