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 441, Wed Jan 19 11:22:07 2005 UTC revision 442, Wed Jan 19 11:25:39 2005 UTC
# Line 133  Line 133 
133  }  }
134    
135  /**  /**
136   * Replace the structure of C by the structure of CA   * Replace the structure of C by the structure of CL^{-1} where L is the
137     * unit lower triangular sparse matrix from an LDL' Cholesky decomposition
138   *   *
139   * @param A a unit lower triangular cscBlocked object   * @param anc number of columns in A
140     * @param Parent parent array for A
141   * @param C a cscBlocked object to be updated   * @param C a cscBlocked object to be updated
142   */   */
143  static void  static void
144  symbolic_right_unit_mm(SEXP A, SEXP C)  symbolic_right_unit_sm(int anc, const int Parent[], SEXP C)
145  {  {
146      SEXP aip = GET_SLOT(A, Matrix_iSym),      SEXP cip = GET_SLOT(C, Matrix_iSym),
         app = GET_SLOT(A, Matrix_pSym),  
         cip = GET_SLOT(C, Matrix_iSym),  
147          cpp = GET_SLOT(C, Matrix_pSym);          cpp = GET_SLOT(C, Matrix_pSym);
148      int *Flag,      int *Flag,
         *ai = INTEGER(aip),  
         *ap = INTEGER(app),  
149          *ci = INTEGER(cip),          *ci = INTEGER(cip),
150          *cp = INTEGER(cpp),          *cp = INTEGER(cpp),
151          *ncp,          *ncp,
         anc = length(app) - 1,  
         anz = length(aip),  
152          cnr, cnz = length(cip),          cnr, cnz = length(cip),
153          i, j;          i, j;
154    
155      if ((length(cpp) - 1) != anc) /* A is square so can compare no of cols */      if ((length(cpp) - 1) != anc) /* A is square so can compare no of cols */
156          error("No. of rows in A (%d) does not match no. of cols in C (%d)",          error("No. of rows in A (%d) does not match no. of cols in C (%d)",
157                anc, length(cpp) - 1);                anc, length(cpp) - 1);
158      if (anz < 1) return;        /* A is the identity */      i = 1;                      /* check for A being the identity */
159      cnr = -1;                   /* number of rows in C is max(ci + 1) */      for (j = 0; j < anc; j++) {
160            if (Parent[j] >= 0) {
161                i = 0;
162                break;
163            }
164        }
165        if (i) return;              /* A is the identity */
166        cnr = 0;                    /* number of rows in C (= max(ci + 1)) */
167      for (i = 0; i < cnz; i++) {      for (i = 0; i < cnz; i++) {
168          int ri = ci[i] + 1;          int ri = ci[i] + 1;
169          if (cnr < ri) cnr = ri;          if (cnr < ri) cnr = ri;
# Line 170  Line 173 
173    
174      ncp[0] = 0;      ncp[0] = 0;
175      for (j = 0; j < anc; j++) {      for (j = 0; j < anc; j++) {
176          int aj2 = ap[j + 1], cj2 = cp[j + 1], ka, kc;          int cj2 = cp[j + 1], kk, kc;
177          for (i = 0; i < cnr; i++) Flag[i] = 0;          for (i = 0; i < cnr; i++) Flag[i] = 0;
178          ncp[j+1] = ncp[j] + cj2 - cp[j];          ncp[j+1] = ncp[j] + cj2 - cp[j];
179                                  /* positions of current column j of C */                                  /* positions of current column j of C */
180          for (kc = cp[j]; kc < cj2; kc++) Flag[ci[kc]] = 1;          for (kc = cp[j]; kc < cj2; kc++) Flag[ci[kc]] = 1;
181                                  /* other positions in column j of product */                                  /* other positions in column j of product */
182          for (ka = ap[j]; ka < aj2; ka++) {          for (kk = Parent[j]; kk >= 0; kk = Parent[kk]) { /* loop over the Parent */
183              int kk = ai[ka], kk2 = cp[kk + 1];              int kk2 = cp[kk + 1];
184              for (kc = cp[kk]; kc < kk2; kc++) {              for (kc = cp[kk]; kc < kk2; kc++) {
185                  if (!Flag[ci[kc]]) {                  if (!Flag[ci[kc]]) {
186                      ncp[j+1]++;                      ncp[j+1]++;
# Line 198  Line 201 
201          for (i = 0; i < nnz; i++) ncx[i] = 1.;          for (i = 0; i < nnz; i++) ncx[i] = 1.;
202                                  /* As Diana Krall said, "Just do it again." */                                  /* As Diana Krall said, "Just do it again." */
203          for (j = 0; j < anc; j++) {          for (j = 0; j < anc; j++) {
204              int aj2 = ap[j + 1], cj2 = cp[j + 1], ka, kc;              int cj2 = cp[j + 1], kc, kk;
205              for (i = 0; i < cnr; i++) Flag[i] = 0;              for (i = 0; i < cnr; i++) Flag[i] = 0;
206              for (kc = cp[j]; kc < cj2; kc++) Flag[ci[kc]] = 1;              for (kc = cp[j]; kc < cj2; kc++) Flag[ci[kc]] = 1;
207              for (ka = ap[j]; ka < aj2; ka++) {              for (kk = Parent[j]; kk >= 0; kk = Parent[kk]) {
208                  int kk = ai[ka], kk2 = cp[kk + 1];                  int kk2 = cp[kk + 1];
209                  for (kc = cp[kk]; kc < kk2; kc++) Flag[ci[kc]] = 1;                  for (kc = cp[kk]; kc < kk2; kc++) Flag[ci[kc]] = 1;
210              }              }
211              for (i = 0; i < cnr; i++) if (Flag[i]) nci[pos++] = i;              for (i = 0; i < cnr; i++) if (Flag[i]) nci[pos++] = i;
# Line 213  Line 216 
216  }  }
217    
218  /**  /**
219   * Replace the structure of C by the structure of CA'   * Replace the structure of C by the structure of CA^{-T}
220   *   *
221   * @param A a unit lower triangular cscBlocked object   * @param A a unit lower triangular cscBlocked object
222   * @param C a cscBlocked object to be updated   * @param C a cscBlocked object to be updated
223   */   */
224  static void  static void
225  symbolic_right_unit_mm_trans(SEXP A, SEXP C)  symbolic_right_unit_mm_trans(int anc, const int Parent[], SEXP C)
226  {  {
227      SEXP aip = GET_SLOT(A, Matrix_iSym),      SEXP cip = GET_SLOT(C, Matrix_iSym),
         app = GET_SLOT(A, Matrix_pSym),  
         cip = GET_SLOT(C, Matrix_iSym),  
228          cpp = GET_SLOT(C, Matrix_pSym);          cpp = GET_SLOT(C, Matrix_pSym);
229      int *ai = INTEGER(aip),      int *ci = INTEGER(cip),
         *ap = INTEGER(app),  
         *ci = INTEGER(cip),  
230          *cp = INTEGER(cpp),          *cp = INTEGER(cpp),
         anc = length(app) - 1,  
         anz = length(aip),  
231          cnz = length(cip),          cnz = length(cip),
232          j, nextra = 0;          i, j, nextra = 0;
233    
234      if ((length(cpp) - 1) != anc)      if ((length(cpp) - 1) != anc)
235          error("No. of cols in A (%d) does not match no. of cols in C (%d)",          error("No. of cols in A (%d) does not match no. of cols in C (%d)",
236                anc, length(cpp) - 1);                anc, length(cpp) - 1);
237      if (anz < 1) return;        /* A is the identity */  
238        i = 1;                      /* check for A being the identity */
239        for (j = 0; j < anc; j++) {
240            if (Parent[j] >= 0) {
241                i = 0;
242                break;
243            }
244        }
245        if (i) return;              /* A is the identity */
246    
247      for (j = 0; j < anc; j++) { /* bound the number of extra triplets */      for (j = 0; j < anc; j++) { /* bound the number of extra triplets */
248          int aj2 = ap[j + 1], cj2 = cp[j + 1], ka, kc;          int cj2 = cp[j + 1], ka, kc;
249          for (ka = ap[j]; ka < aj2; ka++) {          for (ka = Parent[j]; ka >= 0; ka = Parent[ka]) {
250              for (kc = cp[j]; kc < cj2; kc++) {              for (kc = cp[j]; kc < cj2; kc++) {
251                  if (check_csc_index(cp, ci, ai[ka], ci[kc]) < 0) nextra++;                  if (check_csc_index(cp, ci, ci[kc], ka) < 0) nextra++;
252              }              }
253          }          }
254      }      }
# Line 254  Line 260 
260              *Ci = Calloc(ntot, int);              *Ci = Calloc(ntot, int);
261    
262          for (j = 0, pos = cnz; j < anc; j++) {          for (j = 0, pos = cnz; j < anc; j++) {
263              int aj2 = ap[j + 1], cj2 = cp[j + 1], ka, kc;              int cj2 = cp[j + 1], ka, kc;
264              for (ka = ap[j]; ka < aj2; ka++) {              for (ka = Parent[j]; ka >= 0; ka = Parent[ka]) {
265                  for (kc = cp[j]; kc < cj2; kc++) {                  for (kc = cp[j]; kc < cj2; kc++) {
266                      if (check_csc_index(cp, ci, ci[kc], ai[ka]) < 0) {                      if (check_csc_index(cp, ci, ci[kc], ka) < 0) {
267                          Tj[pos] = ai[ka];                          Tj[pos] = ka;
268                          Ti[pos] = ci[kc];                          Ti[pos] = ci[kc];
269                          pos++;                          pos++;
270                      }                      }
271                  }                  }
272              }              }
273          }          }
274          for (j = 0, cnr = -1; j < cnz; j++) {int rr = ci[j]; if (rr > cnr) cnr = rr;}          for (j = 0, cnr = 0; j < cnz; j++) { /* determine number of rows in C */
275          cnr++;                  /* maximum index is 1 less than number of rows */              int rr = ci[j] + 1;
276                if (rr > cnr) cnr = rr;
277            }
278          triplet_to_col(cnr, anc, ntot, Ti, Tj, (double *) NULL,          triplet_to_col(cnr, anc, ntot, Ti, Tj, (double *) NULL,
279                         INTEGER(cpp), Ci, (double *) NULL);                         INTEGER(cpp), Ci, (double *) NULL);
280          cnz = cp[anc];          cnz = cp[anc];
# Line 439  Line 447 
447  void  void
448  lmer_populate(SEXP val)  lmer_populate(SEXP val)
449  {  {
450      SEXP D, L, Linv, Parent, cscB = MAKE_CLASS("cscBlocked"), ZZpO,      SEXP D, L, Parent, cscB = MAKE_CLASS("cscBlocked"), ZZpO,
451          flist = GET_SLOT(val, Matrix_flistSym), perm, Omega,          flist = GET_SLOT(val, Matrix_flistSym), perm, Omega,
452          ZtZ = GET_SLOT(val, Matrix_ZtZSym);          ZtZ = GET_SLOT(val, Matrix_ZtZSym);
453      int j, nf = length(flist);      int j, nf = length(flist);
454      int *nc = INTEGER(GET_SLOT(val, Matrix_ncSym)), *Gp,      int *nc = INTEGER(GET_SLOT(val, Matrix_ncSym)), *Gp,
455          *nlev = Calloc(nf, int), npairs = (nf * (nf + 1))/2;          *nlev = Calloc(nf, int), npairs = (nf * (nf + 1))/2;
456      char *statnms[] = {"factored", "inverted", ""},      char *statnms[] = {"factored", "inverted", ""},
457          *devnms[] = {"ML", "REML", ""};          *devnms[] = {"ML", "REML", ""},
458            *pnms[] = {"index", "block", ""};
459    
460      /* Allocate fixed-sized slots */      /* Allocate fixed-sized slots */
461      SET_SLOT(val, Matrix_statusSym, make_named(LGLSXP, statnms ));      SET_SLOT(val, Matrix_statusSym, make_named(LGLSXP, statnms ));
# Line 459  Line 468 
468      SET_SLOT(val, Matrix_DSym, allocVector(VECSXP, nf));      SET_SLOT(val, Matrix_DSym, allocVector(VECSXP, nf));
469      D = GET_SLOT(val, Matrix_DSym);      D = GET_SLOT(val, Matrix_DSym);
470      setAttrib(D, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));      setAttrib(D, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
     SET_SLOT(val, Matrix_LinvSym, allocVector(VECSXP, nf));  
     Linv = GET_SLOT(val, Matrix_LinvSym);  
     setAttrib(Linv, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));  
471      SET_SLOT(val, Matrix_permSym, allocVector(VECSXP, nf));      SET_SLOT(val, Matrix_permSym, allocVector(VECSXP, nf));
472      perm = GET_SLOT(val, Matrix_permSym);      perm = GET_SLOT(val, Matrix_permSym);
473      setAttrib(perm, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));      setAttrib(perm, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
# Line 493  Line 499 
499      for (j = 0; j < nf; j++) {      for (j = 0; j < nf; j++) {
500          int dind = Lind(j, j), i, k;          int dind = Lind(j, j), i, k;
501          SEXP ctd = VECTOR_ELT(ZZpO, dind); /* diagonal in crosstab */          SEXP ctd = VECTOR_ELT(ZZpO, dind); /* diagonal in crosstab */
502          SEXP Linvj, Ljj, cpp = GET_SLOT(ctd, Matrix_pSym),          SEXP Ljj, cpp = GET_SLOT(ctd, Matrix_pSym),
503              cip = GET_SLOT(ctd, Matrix_iSym);              cip = GET_SLOT(ctd, Matrix_iSym), parent;
504          int *Lp, *Linvp, *Perm, *cp = INTEGER(cpp),          int *Lp, *Perm, *cp = INTEGER(cpp),
505              *ci = INTEGER(cip), *dims, *parent,              *ci = INTEGER(cip), *dims,
506              ncj = length(cpp) - 1,              ncj = length(cpp) - 1,
507              nnz = length(cip);              nnz = length(cip);
         double *dtmp;  
508    
509          SET_VECTOR_ELT(Parent, j, allocVector(INTSXP, ncj));          SET_VECTOR_ELT(Parent, j, make_named(VECSXP, pnms));
510          parent = INTEGER(VECTOR_ELT(Parent, j));          parent = VECTOR_ELT(Parent, j);
511            SET_VECTOR_ELT(parent, 0, allocVector(INTSXP, ncj));
512            SET_VECTOR_ELT(parent, 1, allocVector(INTSXP, ncj));
513          SET_VECTOR_ELT(perm, j, allocVector(INTSXP, ncj));          SET_VECTOR_ELT(perm, j, allocVector(INTSXP, ncj));
514          Perm = INTEGER(VECTOR_ELT(perm, j));          Perm = INTEGER(VECTOR_ELT(perm, j));
515          SET_VECTOR_ELT(L, dind, NEW_OBJECT(cscB));          SET_VECTOR_ELT(L, dind, NEW_OBJECT(cscB));
516          Ljj = VECTOR_ELT(L, dind);          Ljj = VECTOR_ELT(L, dind);
         SET_VECTOR_ELT(Linv, j, NEW_OBJECT(cscB));  
         Linvj = VECTOR_ELT(Linv, j);  
517          SET_SLOT(Ljj, Matrix_pSym, allocVector(INTSXP, ncj + 1));          SET_SLOT(Ljj, Matrix_pSym, allocVector(INTSXP, ncj + 1));
         SET_SLOT(Linvj, Matrix_pSym, allocVector(INTSXP, ncj + 1));  
518          Lp = INTEGER(GET_SLOT(Ljj, Matrix_pSym));          Lp = INTEGER(GET_SLOT(Ljj, Matrix_pSym));
         Linvp = INTEGER(GET_SLOT(Linvj, 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 526  Line 529 
529                                  /* apply to the factor */                                  /* apply to the factor */
530              factor_levels_permute(fac, fcp, Perm, iPerm);              factor_levels_permute(fac, fcp, Perm, iPerm);
531                                  /* symbolic analysis to get Parent */                                  /* symbolic analysis to get Parent */
532              R_ldl_symbolic(ncj, cp, ci, Lp, parent,              R_ldl_symbolic(ncj, cp, ci, Lp, INTEGER(VECTOR_ELT(parent, 0)),
533                           (int *) NULL, (int *) NULL);                           (int *) NULL, (int *) NULL);
534                                  /* decompose the identity to get the row pointers */              for (i = 0; i < ncj; i++)
535              dtmp = REAL(GET_SLOT(ctd, Matrix_xSym));                  INTEGER(VECTOR_ELT(parent, 1))[i] =
536              for (i = 0; i < nnz; i++) dtmp[i] = 0.; /* initialize */                      (INTEGER(VECTOR_ELT(parent, 0))[i] < 0) ? -1 : j;
                                 /* diagonal el is last element in the column */  
             for (i = 0; i < ncj; i++) dtmp[cp[i+1] - 1] = 1.;  
537              nnz = Lp[ncj];              nnz = Lp[ncj];
538              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, nnz));              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, nnz));
539              SET_SLOT(Ljj, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], nnz));              SET_SLOT(Ljj, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], nnz));
             i = R_ldl_numeric(ncj, cp, ci, dtmp, Lp, parent,  
                                INTEGER(GET_SLOT(Ljj, Matrix_iSym)),  
                                REAL(GET_SLOT(Ljj, Matrix_xSym)),  
                                (double *) R_alloc(ncj, sizeof(double)), /* D */  
                                (int *) NULL, (int *) NULL);  
             if (i < ncj) error("identity matrix is not positive definite?");  
540              Free(iPerm); UNPROTECT(1);              Free(iPerm); UNPROTECT(1);
541          } else {          } else {
542              for (i = 0; i <= ncj; i++) {              for (i = 0; i < ncj; i++) {
543                  Lp[i] = 0;                  Lp[i] = 0;
544                  parent[i] = -1;                  INTEGER(VECTOR_ELT(parent,0))[i] = -1;
545                    INTEGER(VECTOR_ELT(parent,1))[i] = -1;
546                  Perm[i] = i;                  Perm[i] = i;
547              }              }
548                Lp[ncj] = 0;
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          }          }
552                                  /* Derive the diagonal block of Linv */          /* FIXME: Remove the references to Linv and use L instead of ZZpO */
         nnz = parent_inv_ap(ncj, 0, parent, Linvp);  
         SET_SLOT(Linvj, Matrix_iSym, allocVector(INTSXP, nnz));  
         parent_inv_ai(ncj, 0, parent, INTEGER(GET_SLOT(Linvj, Matrix_iSym)));  
         SET_SLOT(Linvj, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], nnz));  
         dtmp = REAL(GET_SLOT(Linvj, Matrix_xSym));  
         for (i = 0; i < nnz; i++) dtmp[i] = 1.;  
553          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 */
554              symbolic_right_unit_mm_trans(Linvj, VECTOR_ELT(ZZpO, Lind(k,j)));              symbolic_right_unit_mm_trans(ncj, INTEGER(VECTOR_ELT(parent, 0)),
555                                             VECTOR_ELT(ZZpO, Lind(k,j)));
556          }          }
557          for (k = j+1; k < nf; k++) { /* Update remaining columns */          for (k = j+1; k < nf; k++) { /* Update remaining columns */
558              for (i = k; i < nf; i++) block_update(ZZpO, j, k, i);              for (i = k; i < nf; i++) block_update(ZZpO, j, k, i);
# Line 569  Line 561 
561              SET_VECTOR_ELT(L, Lind(j,i), duplicate(VECTOR_ELT(ZZpO, Lind(j,i))));              SET_VECTOR_ELT(L, Lind(j,i), duplicate(VECTOR_ELT(ZZpO, Lind(j,i))));
562          }          }
563      }      }
564        /* Convert blockwise Parent arrays to extended Parent arrays */
565        for (j = 0; j < (nf - 1); j++) { /* Parent[nf] does not need conversion */
566            SEXP Ljp1j = VECTOR_ELT(L, Lind(j + 1, j)),
567                LpP = GET_SLOT(Ljp1j, Matrix_pSym);
568            int *Li = INTEGER(GET_SLOT(Ljp1j, Matrix_iSym)),
569                *Lp = INTEGER(LpP),
570                *block = INTEGER(VECTOR_ELT(VECTOR_ELT(Parent, j), 1)),
571                *parent = INTEGER(VECTOR_ELT(VECTOR_ELT(Parent, j), 0)),
572                i, nlev = length(LpP) - 1;
573            for (i = 0; i < nlev; i++) {
574                if (block[i] < 0) {
575                    block[i] = j + 1;
576                    parent[i] = Li[Lp[i]];
577                }
578            }
579        }
580      Free(nlev);      Free(nlev);
581  }  }

Legend:
Removed from v.441  
changed lines
  Added in v.442

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