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 493, Thu Feb 3 18:47:20 2005 UTC revision 536, Thu Feb 10 08:15:12 2005 UTC
# Line 1  Line 1 
1  #include "bCrosstab.h"  #include "bCrosstab.h"
2  /* TODO:  /* TODO:
  * - Use the off-diagonal blocks of L  
  * - Remove the off-diagonal blocks of ZZpO  
3   * - Only do a fill-reducing permutation on the first non-nested factor   * - Only do a fill-reducing permutation on the first non-nested factor
4     * - Alternatively: change the algorithm for the fill-reducing
5     *   permutation to a greedy or a picky algorithm.
6   */   */
7    
8  /**  /**
  * Allocate a 3-dimensional array  
  *  
  * @param TYP The R Type code (e.g. INTSXP)  
  * @param nr number of rows  
  * @param nc number of columns  
  * @param nf number of faces  
  *  
  * @return A 3-dimensional array of the indicated dimensions and type  
  */  
 static  
 SEXP alloc3Darray(int TYP, int nr, int nc, int nf)  
 {  
     SEXP val, dd = PROTECT(allocVector(INTSXP, 3));  
   
     INTEGER(dd)[0] = nr; INTEGER(dd)[1] = nc; INTEGER(dd)[2] = nf;  
     val = allocArray(TYP, dd);  
     UNPROTECT(1);  
     return val;  
 }  
   
 /**  
  * Calculate the zero-based index in a row-wise packed lower triangular matrix.  
  * This is used for the arrays of blocked sparse matrices.  
  *  
  * @param i column number (zero-based)  
  * @param k row number (zero-based)  
  *  
  * @return The index of the (k,i) element of a packed lower triangular matrix  
  */  
 static R_INLINE  
 int Lind(int k, int i)  
 {  
     if (k < i) error("Lind(k = %d, i = %d) must have k >= i", k, i);  
     return (k * (k + 1))/2 + i;  
 }  
   
 /**  
9   * Apply a permutation to an index vector   * Apply a permutation to an index vector
10   *   *
11   * @param i vector of 0-based indices   * @param i vector of 0-based indices
# Line 77  Line 40 
40  }  }
41    
42  /**  /**
  * Check for the existence of the (row, col) pair in a csc structure.  
  *  
  * @param p vector of column pointers  
  * @param i vector of row indices  
  * @param row row index  
  * @param col column index  
  *  
  * @return index of element at (row, col) if it exists, otherwise -1  
  */  
 static R_INLINE int  
 check_csc_index(const int p[], const int i[], int row, int col)  
 {  
     int k, k2 = p[col + 1];  
                                 /* use a linear search for now */  
                                 /* perhaps replace by bsearch later */  
     for (k = p[col]; k < k2; k++) {  
         if (i[k] == row) return k;  
     }  
     return -1;  
 }  
   
 static int*  
 expand_column_pointers(int ncol, const int mp[], int mj[])  
 {  
     int j;  
     for (j = 0; j < ncol; j++) {  
         int j2 = mp[j+1], jj;  
         for (jj = mp[j]; jj < j2; jj++) mj[jj] = j;  
     }  
     return mj;  
 }  
   
 /**  
43   * Replace the structure of C by the structure of CL^{-1} where L is the   * Replace the structure of C by the structure of CL^{-1} where L is the
44   * unit lower triangular sparse matrix from an LDL' Cholesky decomposition   * unit lower triangular sparse matrix from an LDL' Cholesky decomposition
45   *   *
# Line 156  Line 86 
86                                  /* positions of current column j of C */                                  /* positions of current column j of C */
87          for (kc = cp[j]; kc < cj2; kc++) Flag[ci[kc]] = 1;          for (kc = cp[j]; kc < cj2; kc++) Flag[ci[kc]] = 1;
88                                  /* other positions in column j of product */                                  /* other positions in column j of product */
89          for (kk = Parent[j]; kk >= 0; kk = Parent[kk]) { /* loop over the Parent */          for (kk = Parent[j]; kk >= 0; kk = Parent[kk]) {
90              int kk2 = cp[kk + 1];              int kk2 = cp[kk + 1];
91              for (kc = cp[kk]; kc < kk2; kc++) {              for (kc = cp[kk]; kc < kk2; kc++) {
92                  if (!Flag[ci[kc]]) {                  if (!Flag[ci[kc]]) {
# Line 225  Line 155 
155          int cj2 = cp[j + 1], ka, kc;          int cj2 = cp[j + 1], ka, kc;
156          for (ka = Parent[j]; ka >= 0; ka = Parent[ka]) {          for (ka = Parent[j]; ka >= 0; ka = Parent[ka]) {
157              for (kc = cp[j]; kc < cj2; kc++) {              for (kc = cp[j]; kc < cj2; kc++) {
158                  if (check_csc_index(cp, ci, ci[kc], ka) < 0) nextra++;                  if (check_csc_index(cp, ci, ci[kc], ka, -1) < 0) nextra++;
159              }              }
160          }          }
161      }      }
# Line 240  Line 170 
170              int cj2 = cp[j + 1], ka, kc;              int cj2 = cp[j + 1], ka, kc;
171              for (ka = Parent[j]; ka >= 0; ka = Parent[ka]) {              for (ka = Parent[j]; ka >= 0; ka = Parent[ka]) {
172                  for (kc = cp[j]; kc < cj2; kc++) {                  for (kc = cp[j]; kc < cj2; kc++) {
173                      if (check_csc_index(cp, ci, ci[kc], ka) < 0) {                      if (check_csc_index(cp, ci, ci[kc], ka, -1) < 0) {
174                          Tj[pos] = ka;                          Tj[pos] = ka;
175                          Ti[pos] = ci[kc];                          Ti[pos] = ci[kc];
176                          pos++;                          pos++;
# Line 296  Line 226 
226          int i1, kk, i2 = ip[jj + 1], k2 = kp[jj + 1];          int i1, kk, i2 = ip[jj + 1], k2 = kp[jj + 1];
227          for (kk = kp[jj]; kk < k2; kk++) {          for (kk = kp[jj]; kk < k2; kk++) {
228              for (i1 = ip[jj]; i1 < i2; i1++) {              for (i1 = ip[jj]; i1 < i2; i1++) {
229                      if ((check_csc_index(tp, ti, ii[i1], ki[kk]) < 0) &&                      if ((check_csc_index(tp, ti, ii[i1], ki[kk], -1) < 0) &&
230                                  /* only update upper triangle of                                  /* only update upper triangle of
231                                   * diagonal blocks */                                   * diagonal blocks */
232                          ((k != i) || (ii[i1] <= ki[kk]))) extra++;                          ((k != i) || (ii[i1] <= ki[kk]))) extra++;
# Line 323  Line 253 
253              int i1, kk, i2 = ip[jj + 1], k2 = kp[jj + 1];              int i1, kk, i2 = ip[jj + 1], k2 = kp[jj + 1];
254              for (kk = kp[jj]; kk < k2; kk++) {              for (kk = kp[jj]; kk < k2; kk++) {
255                  for (i1 = ip[jj]; i1 < i2; i1++) {                  for (i1 = ip[jj]; i1 < i2; i1++) {
256                      if ((check_csc_index(tp, ti, ii[i1], ki[kk]) < 0) &&                      if ((check_csc_index(tp, ti, ii[i1], ki[kk], -1) < 0) &&
257                          ((k != i) || (ii[i1] <= ki[kk]))) {                          ((k != i) || (ii[i1] <= ki[kk]))) {
258                          Ti[pos] = ii[i1];                          Ti[pos] = ii[i1];
259                          Tj[pos] = ki[kk];                          Tj[pos] = ki[kk];
# Line 332  Line 262 
262                  }                  }
263              }              }
264          }          }
265          /* Pass nlev instead of doing this.  The dimensions are nlev[i], nlev[k] */          /* FIXME: Pass nlev instead -  dimensions are nlev[i], nlev[k] */
266          /* Determine maximum row index in T */          /* Determine maximum row index in T */
267          tnr = -1; for (jj = 0; jj < ntot; jj++) if (Ti[jj] > tnr) tnr = Ti[jj];          tnr = -1; for (jj = 0; jj < ntot; jj++) if (Ti[jj] > tnr) tnr = Ti[jj];
268          tnr++;                  /* increment by 1 to get number of rows */          tnr++;                  /* increment by 1 to get number of rows */
# Line 342  Line 272 
272          SET_SLOT(tb, Matrix_iSym, allocVector(INTSXP, nnz));          SET_SLOT(tb, Matrix_iSym, allocVector(INTSXP, nnz));
273          Memcpy(INTEGER(GET_SLOT(tb, Matrix_iSym)), Ai, nnz);          Memcpy(INTEGER(GET_SLOT(tb, Matrix_iSym)), Ai, nnz);
274          dims = INTEGER(getAttrib(GET_SLOT(tb, Matrix_xSym), R_DimSymbol));          dims = INTEGER(getAttrib(GET_SLOT(tb, Matrix_xSym), R_DimSymbol));
275          SET_SLOT(tb, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], nnz));          SET_SLOT(tb, Matrix_xSym,
276                     alloc3Darray(REALSXP, dims[0], dims[1], nnz));
277          Ax = REAL(GET_SLOT(tb, Matrix_xSym));          Ax = REAL(GET_SLOT(tb, Matrix_xSym));
278          for (j = 0; j < nnz; j++) Ax[j] = 1.;          for (j = 0; j < nnz; j++) Ax[j] = 1.;
279          Free(Ai); Free(Ti); Free(Tj);          Free(Ai); Free(Ti); Free(Tj);
# Line 361  Line 292 
292   * @param pperm inverse of the permutation   * @param pperm inverse of the permutation
293   */   */
294  static void  static void
295  bCrosstab_permute(SEXP ctab, int nf, int jj, const int nlev[], const int iperm[])  bCrosstab_permute(SEXP ctab, int nf, int jj,
296                      const int nlev[], const int iperm[])
297  {  {
298      int j;      int j;
299      for (j = 0; j < nf; j++) {      for (j = 0; j < nf; j++) {
# Line 443  Line 375 
375  {  {
376      SEXP D, L, Parent, ZZpO, flist = GET_SLOT(val, Matrix_flistSym),      SEXP D, L, Parent, ZZpO, flist = GET_SLOT(val, Matrix_flistSym),
377          perm, Omega, ZtZ = GET_SLOT(val, Matrix_ZtZSym);          perm, Omega, ZtZ = GET_SLOT(val, Matrix_ZtZSym);
378        SEXP fnms = getAttrib(flist, R_NamesSymbol);
379      int j, k, nf = length(flist);      int j, k, nf = length(flist);
380      int *nc = INTEGER(GET_SLOT(val, Matrix_ncSym)), *Gp,      int *nc = INTEGER(GET_SLOT(val, Matrix_ncSym)), *Gp,
381          *nlev = Calloc(nf, int), npairs = (nf * (nf + 1))/2;          *nlev = Calloc(nf, int), npairs = (nf * (nf + 1))/2;
# Line 455  Line 388 
388      SET_SLOT(val, Matrix_devianceSym, Matrix_make_named(REALSXP, devnms));      SET_SLOT(val, Matrix_devianceSym, Matrix_make_named(REALSXP, devnms));
389      SET_SLOT(val, Matrix_devCompSym, allocVector(REALSXP, 4));      SET_SLOT(val, Matrix_devCompSym, allocVector(REALSXP, 4));
390      /* Allocate slots that are lists of length nf */      /* Allocate slots that are lists of length nf */
391      SET_SLOT(val, Matrix_ZZpOSym, allocVector(VECSXP, nf));      ZZpO = ALLOC_SLOT(val, Matrix_ZZpOSym, VECSXP, nf);
392      ZZpO = GET_SLOT(val, Matrix_ZZpOSym);      setAttrib(ZZpO, R_NamesSymbol, duplicate(fnms));
393      setAttrib(ZZpO, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));      D = ALLOC_SLOT(val, Matrix_DSym, VECSXP, nf);
394      SET_SLOT(val, Matrix_DSym, allocVector(VECSXP, nf));      setAttrib(D, R_NamesSymbol, duplicate(fnms));
395      D = GET_SLOT(val, Matrix_DSym);      perm = ALLOC_SLOT(val, Matrix_permSym, VECSXP, nf);
396      setAttrib(D, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));      setAttrib(perm, R_NamesSymbol, duplicate(fnms));
397      SET_SLOT(val, Matrix_permSym, allocVector(VECSXP, nf));      Parent = ALLOC_SLOT(val, Matrix_ParentSym, VECSXP, nf);
398      perm = GET_SLOT(val, Matrix_permSym);      setAttrib(Parent, R_NamesSymbol, duplicate(fnms));
399      setAttrib(perm, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));      Omega = ALLOC_SLOT(val, Matrix_OmegaSym, VECSXP, nf);
400      SET_SLOT(val, Matrix_ParentSym, allocVector(VECSXP, nf));      setAttrib(Omega, R_NamesSymbol, duplicate(fnms));
     Parent = GET_SLOT(val, Matrix_ParentSym);  
     setAttrib(Parent, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));  
     SET_SLOT(val, Matrix_OmegaSym, allocVector(VECSXP, nf));  
     Omega = GET_SLOT(val, Matrix_OmegaSym);  
     setAttrib(Omega, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));  
401    
402      /* Allocate peculiar length slots */      /* Allocate peculiar length slots */
403      SET_SLOT(val, Matrix_LSym, allocVector(VECSXP, npairs));      SET_SLOT(val, Matrix_LSym, allocVector(VECSXP, npairs));
# Line 484  Line 412 
412          SET_VECTOR_ELT(Omega, j, allocMatrix(REALSXP, nc[j], nc[j]));          SET_VECTOR_ELT(Omega, j, allocMatrix(REALSXP, nc[j], nc[j]));
413          SET_VECTOR_ELT(ZZpO, j, duplicate(VECTOR_ELT(ZtZ, Lind(j, j))));          SET_VECTOR_ELT(ZZpO, j, duplicate(VECTOR_ELT(ZtZ, Lind(j, j))));
414          for (k = j; k < nf; k++)          for (k = j; k < nf; k++)
415              SET_VECTOR_ELT(L, Lind(k, j), duplicate(VECTOR_ELT(ZtZ, Lind(k, j))));              SET_VECTOR_ELT(L, Lind(k, j),
416                               duplicate(VECTOR_ELT(ZtZ, Lind(k, j))));
417      }      }
418      SET_SLOT(val, Matrix_XtXSym, allocMatrix(REALSXP, nc[nf], nc[nf]));      SET_SLOT(val, Matrix_XtXSym, allocMatrix(REALSXP, nc[nf], nc[nf]));
419      AZERO(REAL(GET_SLOT(val, Matrix_XtXSym)), nc[nf] * nc[nf]);      AZERO(REAL(GET_SLOT(val, Matrix_XtXSym)), nc[nf] * nc[nf]);
# Line 531  Line 460 
460                      (INTEGER(VECTOR_ELT(parent, 0))[i] < 0) ? -1 : j;                      (INTEGER(VECTOR_ELT(parent, 0))[i] < 0) ? -1 : j;
461              nnz = Lp[ncj];              nnz = Lp[ncj];
462              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, nnz));              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, nnz));
463              SET_SLOT(Ljj, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], nnz));              SET_SLOT(Ljj, Matrix_xSym,
464                         alloc3Darray(REALSXP, dims[0], dims[1], nnz));
465              Free(iPerm); UNPROTECT(1);              Free(iPerm); UNPROTECT(1);
466          } else {          } else {
467              for (i = 0; i < ncj; i++) {              for (i = 0; i < ncj; i++) {
# Line 542  Line 472 
472              }              }
473              Lp[ncj] = 0;              Lp[ncj] = 0;
474              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, 0));              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, 0));
475              SET_SLOT(Ljj, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], 0));              SET_SLOT(Ljj, Matrix_xSym,
476                         alloc3Darray(REALSXP, dims[0], dims[1], 0));
477          }          }
478          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 */
479              symbolic_right_unit_mm_trans(ncj, INTEGER(VECTOR_ELT(parent, 0)),              symbolic_right_unit_mm_trans(ncj, INTEGER(VECTOR_ELT(parent, 0)),

Legend:
Removed from v.493  
changed lines
  Added in v.536

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