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 443, Wed Jan 19 11:28:00 2005 UTC revision 444, Wed Jan 19 19:05:15 2005 UTC
# Line 447  Line 447 
447  void  void
448  lmer_populate(SEXP val)  lmer_populate(SEXP val)
449  {  {
450      SEXP D, L, Parent, cscB = MAKE_CLASS("cscBlocked"), ZZpO,      SEXP D, L, Parent, ZZpO, flist = GET_SLOT(val, Matrix_flistSym),
451          flist = GET_SLOT(val, Matrix_flistSym), perm, Omega,          perm, Omega, ZtZ = GET_SLOT(val, Matrix_ZtZSym);
452          ZtZ = GET_SLOT(val, Matrix_ZtZSym);      int j, k, nf = length(flist);
     int j, nf = length(flist);  
453      int *nc = INTEGER(GET_SLOT(val, Matrix_ncSym)), *Gp,      int *nc = INTEGER(GET_SLOT(val, Matrix_ncSym)), *Gp,
454          *nlev = Calloc(nf, int), npairs = (nf * (nf + 1))/2;          *nlev = Calloc(nf, int), npairs = (nf * (nf + 1))/2;
455      char *statnms[] = {"factored", "inverted", ""},      char *statnms[] = {"factored", "inverted", ""},
# Line 461  Line 460 
460      SET_SLOT(val, Matrix_statusSym, make_named(LGLSXP, statnms ));      SET_SLOT(val, Matrix_statusSym, make_named(LGLSXP, statnms ));
461      SET_SLOT(val, Matrix_devianceSym, make_named(REALSXP, devnms));      SET_SLOT(val, Matrix_devianceSym, make_named(REALSXP, devnms));
462      SET_SLOT(val, Matrix_devCompSym, allocVector(REALSXP, 4));      SET_SLOT(val, Matrix_devCompSym, allocVector(REALSXP, 4));
     /* Copy ZtZ to ZZpO */  
     SET_SLOT(val, Matrix_ZZpOSym, duplicate(ZtZ));  
     ZZpO = GET_SLOT(val, Matrix_ZZpOSym);  
463      /* Allocate slots that are lists of length nf */      /* Allocate slots that are lists of length nf */
464        SET_SLOT(val, Matrix_ZZpOSym, allocVector(VECSXP, nf));
465        ZZpO = GET_SLOT(val, Matrix_ZZpOSym);
466        setAttrib(ZZpO, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
467      SET_SLOT(val, Matrix_DSym, allocVector(VECSXP, nf));      SET_SLOT(val, Matrix_DSym, allocVector(VECSXP, nf));
468      D = GET_SLOT(val, Matrix_DSym);      D = GET_SLOT(val, Matrix_DSym);
469      setAttrib(D, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));      setAttrib(D, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
# Line 489  Line 488 
488          Gp[j + 1] = Gp[j] + nc[j] * nlev[j];          Gp[j + 1] = Gp[j] + nc[j] * nlev[j];
489          SET_VECTOR_ELT(D, j, alloc3Darray(REALSXP, nc[j], nc[j], nlev[j]));          SET_VECTOR_ELT(D, j, alloc3Darray(REALSXP, nc[j], nc[j], nlev[j]));
490          SET_VECTOR_ELT(Omega, j, allocMatrix(REALSXP, nc[j], nc[j]));          SET_VECTOR_ELT(Omega, j, allocMatrix(REALSXP, nc[j], nc[j]));
491            SET_VECTOR_ELT(ZZpO, j, duplicate(VECTOR_ELT(ZtZ, Lind(j, j))));
492            for (k = j; k < nf; k++)
493                SET_VECTOR_ELT(L, Lind(k, j), duplicate(VECTOR_ELT(ZtZ, Lind(k, j))));
494      }      }
495      SET_SLOT(val, Matrix_XtXSym, allocMatrix(REALSXP, nc[nf], nc[nf]));      SET_SLOT(val, Matrix_XtXSym, allocMatrix(REALSXP, nc[nf], nc[nf]));
496      AZERO(REAL(GET_SLOT(val, Matrix_XtXSym)), nc[nf] * nc[nf]);      AZERO(REAL(GET_SLOT(val, Matrix_XtXSym)), nc[nf] * nc[nf]);
# Line 497  Line 499 
499      SET_SLOT(val, Matrix_ZtXSym, allocMatrix(REALSXP, Gp[nf], nc[nf]));      SET_SLOT(val, Matrix_ZtXSym, allocMatrix(REALSXP, Gp[nf], nc[nf]));
500      SET_SLOT(val, Matrix_RZXSym, allocMatrix(REALSXP, Gp[nf], nc[nf]));      SET_SLOT(val, Matrix_RZXSym, allocMatrix(REALSXP, Gp[nf], nc[nf]));
501      for (j = 0; j < nf; j++) {      for (j = 0; j < nf; j++) {
502          int dind = Lind(j, j), i, k;          int dind = Lind(j, j), i;
503          SEXP ctd = VECTOR_ELT(ZZpO, dind); /* diagonal in crosstab */          SEXP ctd = VECTOR_ELT(ZZpO, j); /* diagonal in crosstab */
504          SEXP Ljj, cpp = GET_SLOT(ctd, Matrix_pSym),          SEXP Ljj = VECTOR_ELT(L, dind),
505                cpp = GET_SLOT(ctd, Matrix_pSym),
506              cip = GET_SLOT(ctd, Matrix_iSym), parent;              cip = GET_SLOT(ctd, Matrix_iSym), parent;
507          int *Lp, *Perm, *cp = INTEGER(cpp),          int *Lp = INTEGER(GET_SLOT(Ljj, Matrix_pSym)), *Perm,
508                *cp = INTEGER(cpp),
509              *ci = INTEGER(cip), *dims,              *ci = INTEGER(cip), *dims,
510              ncj = length(cpp) - 1,              ncj = length(cpp) - 1,
511              nnz = length(cip);              nnz = length(cip);
# Line 512  Line 516 
516          SET_VECTOR_ELT(parent, 1, allocVector(INTSXP, ncj));          SET_VECTOR_ELT(parent, 1, allocVector(INTSXP, ncj));
517          SET_VECTOR_ELT(perm, j, allocVector(INTSXP, ncj));          SET_VECTOR_ELT(perm, j, allocVector(INTSXP, ncj));
518          Perm = INTEGER(VECTOR_ELT(perm, j));          Perm = INTEGER(VECTOR_ELT(perm, j));
         SET_VECTOR_ELT(L, dind, NEW_OBJECT(cscB));  
         Ljj = VECTOR_ELT(L, dind);  
         SET_SLOT(Ljj, Matrix_pSym, allocVector(INTSXP, ncj + 1));  
         Lp = INTEGER(GET_SLOT(Ljj, Matrix_pSym));  
519          dims = INTEGER(getAttrib(GET_SLOT(ctd, Matrix_xSym), R_DimSymbol));          dims = INTEGER(getAttrib(GET_SLOT(ctd, Matrix_xSym), R_DimSymbol));
520          if (nnz > ncj) {        /* calculate fill-reducing permutation */          if (nnz > ncj) {        /* calculate fill-reducing permutation */
521              SEXP fac = VECTOR_ELT(flist, j);              SEXP fac = VECTOR_ELT(flist, j);
# Line 549  Line 549 
549              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, 0));              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, 0));
550              SET_SLOT(Ljj, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], 0));              SET_SLOT(Ljj, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], 0));
551          }          }
         /* FIXME: Remove the references to Linv and use L instead of ZZpO */  
552          for (k = j+1; k < nf; k++) { /* Update other blocks in this column */          for (k = j+1; k < nf; k++) { /* Update other blocks in this column */
553              symbolic_right_unit_mm_trans(ncj, INTEGER(VECTOR_ELT(parent, 0)),              symbolic_right_unit_mm_trans(ncj, INTEGER(VECTOR_ELT(parent, 0)),
554                                           VECTOR_ELT(ZZpO, Lind(k,j)));                                           VECTOR_ELT(L, Lind(k,j)));
555          }          }
556          for (k = j+1; k < nf; k++) { /* Update remaining columns */          for (k = j+1; k < nf; k++) { /* Update remaining columns */
557              for (i = k; i < nf; i++) block_update(ZZpO, j, k, i);              for (i = k; i < nf; i++) block_update(L, j, k, i);
         }  
         for (i= 0; i < j; i++) { /* copy blocks to the left */  
             SET_VECTOR_ELT(L, Lind(j,i), duplicate(VECTOR_ELT(ZZpO, Lind(j,i))));  
558          }          }
559      }      }
560      /* Convert blockwise Parent arrays to extended Parent arrays */      /* Convert blockwise Parent arrays to extended Parent arrays */

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

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