# SCM Repository

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

# Diff of /pkg/src/bCrosstab.c

revision 359, Sat Nov 27 17:58:23 2004 UTC revision 370, Sun Dec 5 13:52:06 2004 UTC
# Line 10  Line 10
10   * @return The 0-based index of the (i,k) element of a row-wise packed lower   * @return The 0-based index of the (i,k) element of a row-wise packed lower
11   * triangular matrix.   * triangular matrix.
12   */   */
13  static R_INLINE  static R_INLINE int
14  int Lind(int i, int k)  Lind(int i, int k)
15  {  {
16      return (i * (i + 1))/2 + k;      return (i * (i + 1))/2 + k;
17  }  }
18
19  /**  /**
20   * Permute an index vector   * Apply a permutation to an index vector
21   *   *
22   * @param i vector of 0-based indices   * @param i vector of 0-based indices
23   * @param nnz length of vector i   * @param nnz length of vector i
24   * @param perm 0-based permutation vector of length max(i)   * @param perm 0-based permutation vector of length max(i)
25   */   */
26  static R_INLINE  static R_INLINE void
27  void ind_permute(int i[], int nnz, const int perm[])  ind_permute(int i[], int nnz, const int perm[])
28  {  {
29      int j;      int j;
30      for (j = 0; j < nnz; j++) i[j] = perm[i[j]];      for (j = 0; j < nnz; j++) i[j] = perm[i[j]];
# Line 37  Line 37
37   * @param j vector of 0-based column indices   * @param j vector of 0-based column indices
38   * @param nnz length of index vectors   * @param nnz length of index vectors
39   */   */
40  static R_INLINE  static R_INLINE void
41  void make_upper_triangular(int i[], int j[], int nnz)  make_upper_triangular(int i[], int j[], int nnz)
42  {  {
43      int k;      int k;
44      for (k = 0; k < nnz; k++) {      for (k = 0; k < nnz; k++) {
# Line 58  Line 58
58   *   *
59   * @return pointer to a named list of length n   * @return pointer to a named list of length n
60   */   */
61  static  static SEXP
62  SEXP make_named_list(int n, char **names)  make_named_list(int n, char **names)
63  {  {
64      SEXP ans = PROTECT(allocVector(VECSXP, n)),      SEXP ans = PROTECT(allocVector(VECSXP, n)),
65          nms = PROTECT(allocVector(STRSXP, n));          nms = PROTECT(allocVector(STRSXP, n));
# Line 71  Line 71
71      return ans;      return ans;
72  }  }
73
74  static R_INLINE  /**
75  int check_csc_index(const int p[], const int i[], int row, int col)   * Check for the existence of the (row, col) pair in a csc structure.
76     *
77     * @param p vector of column pointers
78     * @param i vector of row indices
79     * @param row row index
80     * @param col column index
81     *
82     * @return index of element at (row, col) if it exists, otherwise -1
83     */
84    static R_INLINE int
85    check_csc_index(const int p[], const int i[], int row, int col)
86  {  {
87      int k, k2 = p[col + 1];      int k, k2 = p[col + 1];
88                                  /* use a linear search for now */                                  /* use a linear search for now */
# Line 84  Line 94
94  }  }
95
96  /**  /**
97   * Update a diagonal block   * Replace the structure of C by the structure of CA
98   *   *
99   * @param ctab pointer to a blocked crosstabulation object   * @param A a unit lower tscMatrix
100   * @param j index of updating column   * @param C cscMatrix to be updated
* @param i index of diagonal block to be updated
101   */   */
102  static void diag_update(SEXP ctab, int j, int i)  static void
103    symbolic_right_unit_mm(SEXP A, SEXP C)
104  {  {
105      SEXP db = VECTOR_ELT(ctab, Lind(i, i)),      SEXP aip = GET_SLOT(A, Matrix_iSym),
106          jb = VECTOR_ELT(ctab, Lind(i, j));          app = GET_SLOT(A, Matrix_pSym),
107      SEXP dpp = GET_SLOT(db, Matrix_pSym),          cip = GET_SLOT(C, Matrix_iSym),
108          jpp = GET_SLOT(jb, Matrix_pSym);          cpp = GET_SLOT(C, Matrix_pSym);
109      int *di = INTEGER(GET_SLOT(db, Matrix_iSym)),      int *Flag,
110          *dp = INTEGER(dpp),          *ai = INTEGER(aip),
111          *ji = INTEGER(GET_SLOT(jb, Matrix_iSym)),          *ap = INTEGER(app),
112          *jp = INTEGER(jpp),          *ci = INTEGER(cip),
113          dnc = length(dpp) - 1,          *cp = INTEGER(cpp),
114          jnc = length(jpp) - 1;          *ncp,
115      int jj, extra;          anc = length(app) - 1,
116                                  /* bound the number of extra elements */          anz = length(aip),
117      extra = 0;          cnr, cnz = length(cip),
118      for (jj = 0; jj < jnc; jj++) {          i, j;
119          int k, kk, k2 = jp[jj + 1];
120          for (k = jp[jj]; k < k2; k++) {      if ((length(cpp) - 1) != anc)
121              for (kk = k; kk < k2; kk++) {          error("No. of rows in A (%d) does not match no. of cols in C (%d)",
122                  if (check_csc_index(dp, di, ji[k], ji[kk]) < 0) extra++;                anc, length(cpp) - 1);
123        if (anz < 1) return;        /* C is the identity */
124        cnr = -1;                   /* number of rows in C is max(ci + 1) */
125        for (i = 0; i < cnz; i++) {
126            int ri = ci[i] + 1;
127            if (cnr < ri) cnr = ri;
128        }
129        Flag = Calloc(cnr, int);
130        ncp = Calloc(anc + 1, int); /* new column pointers */
131
132        ncp[0] = 0;
133        for (j = 0; j < anc; j++) {
134            int aj2 = ap[j + 1], cj2 = cp[j + 1], ka, kc;
135            for (i = 0; i < cnr; i++) Flag[i] = 0;
136            ncp[j+1] = ncp[j] + cj2 - cp[j];
137                                    /* positions of current column j of C */
138            for (kc = cp[j]; kc < cj2; kc++) Flag[ci[kc]] = 1;
139                                    /* other positions in column j of product */
140            for (ka = ap[j]; ka < aj2; ka++) {
141                int kk = ai[ka], kk2 = cp[kk + 1];
142                for (kc = cp[kk]; kc < kk2; kc++) {
143                    if (!Flag[ci[kc]]) {
144                        ncp[j+1]++;
145                        Flag[ci[kc]] = 1;
146                    }
147                }
148            }
149        }
150        if (ncp[anc] > cp[anc]) {
151            int *nci, nnz = ncp[anc], pos = 0;
152            double *ncx;
153
154            SET_SLOT(C, Matrix_iSym, allocVector(INTSXP,  nnz));
155            nci = INTEGER(GET_SLOT(C, Matrix_iSym));
156            SET_SLOT(C, Matrix_xSym, allocVector(REALSXP,  nnz));
157            ncx = REAL(GET_SLOT(C, Matrix_xSym));
158            for (i = 0; i < nnz; i++) ncx[i] = 1.;
159                                    /* As Diana Krall said, "Just do it again." */
160            for (j = 0; j < anc; j++) {
161                int aj2 = ap[j + 1], cj2 = cp[j + 1], ka, kc;
162                for (i = 0; i < cnr; i++) Flag[i] = 0;
163                for (kc = cp[j]; kc < cj2; kc++) Flag[ci[kc]] = 1;
164                for (ka = ap[j]; ka < aj2; ka++) {
165                    int kk = ai[ka], kk2 = cp[kk + 1];
166                    for (kc = cp[kk]; kc < kk2; kc++) Flag[ci[kc]] = 1;
167              }              }
168                for (i = 0; i < cnr; i++) if (Flag[i]) nci[pos++] = i;
169          }          }
170            Memcpy(cp, ncp, anc + 1);
171      }      }
172      if (!extra) return;      Free(Flag); Free(ncp);
{
int pos, nnz = dp[dnc];
int ntot = nnz + extra;
int *Ai = Calloc(ntot, int),
*Ti = Calloc(ntot, int),
*Tj = Calloc(ntot, int);
double *Ax;

Memcpy(Ti, di, nnz);    /* make a copy of the row indices */
pos = 0;                /* fill in the column indices */
for (jj = 0; jj < dnc; jj++) {
int j2 = dp[jj + 1];
while (pos < j2) {
Tj[pos] = jj;
pos++;
}
}
/* add the extra elements */
for (jj = 0; jj < jnc; jj++) {
int k, kk, k2 = jp[jj + 1];
for (k = jp[jj]; k < k2; k++) {
for (kk = k; kk < k2; kk++) {
if (check_csc_index(dp, di, ji[k], ji[kk]) < 0) {
Ti[pos] = ji[k];
Tj[pos] = ji[kk];
pos++;
}
}
}
}
triplet_to_col(dnc, dnc, ntot, Ti, Tj, (double *) NULL,
dp, Ai, (double *) NULL);
nnz = dp[dnc];
SET_SLOT(db, Matrix_iSym, allocVector(INTSXP, nnz));
Memcpy(INTEGER(GET_SLOT(db, Matrix_iSym)), Ai, nnz);
SET_SLOT(db, Matrix_xSym, allocVector(REALSXP, nnz));
Ax = REAL(GET_SLOT(db, Matrix_xSym));
for (j = 0; j < nnz; j++) Ax[j] = 1.;
Free(Ai); Free(Ti); Free(Tj);
return;
}
173  }  }
174
175  /**  /**
176   * Update a block   * Update a block of L in the blocked crosstabulation
177   *   *
178   * @param ctab pointer to a blocked crosstabulation object   * @param ctab pointer to a blocked crosstabulation object
179   * @param j index of updating column   * @param j index of updating column
180   * @param k column index of block to be updated   * @param k column index of block to be updated
181   * @param i row index of block to be updated (j < k <= i)   * @param i row index of block to be updated (j < k <= i)
182   */   */
183  static void block_update(SEXP ctab, int j, int k, int i)  static void
184    block_update(SEXP ctab, int j, int k, int i)
185  {  {
186      SEXP tb = VECTOR_ELT(ctab, Lind(i, k)),      SEXP tb = VECTOR_ELT(ctab, Lind(i, k)),
187          ib = VECTOR_ELT(ctab, Lind(i, j)),          ib = VECTOR_ELT(ctab, Lind(i, j)),
# Line 250  Line 266
266   * @param perm permutation (0-based) to apply   * @param perm permutation (0-based) to apply
267   * @param pperm inverse of the permutation   * @param pperm inverse of the permutation
268   */   */
269  void bCrosstab_permute(SEXP ctab, int nf, int jj,  static void
270    bCrosstab_permute(SEXP ctab, int nf, int jj,
271                         const int perm[], const int iperm[])                         const int perm[], const int iperm[])
272  {  {
273      SEXP trip, ipt, jpt;      SEXP trip, ipt, jpt;
# Line 278  Line 295
295      }      }
296  }  }
297
298    /**
299     * Apply a permutation vector to the levels of a factor.
300     *
301     * The dest pointer is assumed to point to a copy of the src pointer's
302     * contents.
303     *
304     * @param dest pointer to the destination factor
305     * @param src pointer to the source factor
306     * @param perm permutation vector (0-based)
307     * @param iperm inverse permutation vector (0-based)
308     */
309  static void  static void
310  factor_levels_permute(SEXP dest, SEXP src, const int perm[],  factor_levels_permute(SEXP dest, SEXP src, const int perm[],
311                        const int iperm[])                        const int iperm[])
# Line 297  Line 325
325          d[i] = 1 + iperm[s[i]-1];          d[i] = 1 + iperm[s[i]-1];
326  }  }
327
328  SEXP bCrosstab_convert(SEXP bCtab)  /**
329     * Derive the structure of elements of an lmeRep object from a blocked
330     * crosstabulation.
331     *
332     * @param bCtab Pointer to a list containing a blocked crosstabulation
333     * object and the factors used to generate it.
334     *
335     * @return A list of the factor list, the blocked/blocked sparse unit
336     * lower triangular matrix L, the list of diagonal elements of
337     * L-inverse, the permutations and the Parent arrays.
338     */
339    SEXP
340    bCrosstab_convert(SEXP bCtab)
341  {  {
342      char *anms[] = {"flist", "L", "Linv", "perm", "Parent"};      char *anms[] = {"flist", "L", "Linv", "perm", "Parent"};
343      SEXP flist = VECTOR_ELT(bCtab, 0),      SEXP flist = VECTOR_ELT(bCtab, 0),
344          ctab = PROTECT(duplicate(VECTOR_ELT(bCtab, 1))),          ctab = PROTECT(duplicate(VECTOR_ELT(bCtab, 1))),
345          ans = PROTECT(make_named_list(5, anms));          ans = PROTECT(make_named_list(sizeof(anms)/sizeof(char *), anms));
346      SEXP fcp, perm, L, Linv, Parent;      SEXP fcp, perm, L, Linv, Parent;
347      int ctbl = length(ctab), j, nf = length(flist);      int ctbl = length(ctab), j, nf = length(flist);
348
# Line 351  Line 391
391          Lp = INTEGER(GET_SLOT(Ljj, Matrix_pSym));          Lp = INTEGER(GET_SLOT(Ljj, Matrix_pSym));
392          Linvp = INTEGER(GET_SLOT(Linvj, Matrix_pSym));          Linvp = INTEGER(GET_SLOT(Linvj, Matrix_pSym));
393          if (nnz > ncj) {        /* calculate fill-reducing permutation */          if (nnz > ncj) {        /* calculate fill-reducing permutation */
394              int *Li, *Lnz = Calloc(ncj, int), *tmp = Calloc(ncj, int), info;              int *iPerm = Calloc(ncj, int);
395
396              ssc_metis_order(ncj, cp, ci, Perm, tmp);              ssc_metis_order(ncj, cp, ci, Perm, iPerm);
397                                  /* apply to the crosstabulation */                                  /* apply to the crosstabulation */
398              bCrosstab_permute(ctab, nf, j, Perm, tmp);              bCrosstab_permute(ctab, nf, j, Perm, iPerm);
399                                  /* apply to the factor */                                  /* apply to the factor */
400              factor_levels_permute(VECTOR_ELT(fcp, j),              factor_levels_permute(VECTOR_ELT(fcp, j),
401                                    VECTOR_ELT(flist, j), Perm, tmp);                                    VECTOR_ELT(flist, j), Perm, iPerm);
402                                  /* symbolic analysis to get Parent */                                  /* symbolic analysis to get Parent */
403              ldl_symbolic(ncj, cp, ci, Lp, parent, Lnz, tmp,              R_ldl_symbolic(ncj, cp, ci, Lp, parent,
404                           (int *) NULL, (int *) NULL);                           (int *) NULL, (int *) NULL);
405                                  /* decompose the identity to get the row pointers */                                  /* decompose the identity to get the row pointers */
406              dtmp = REAL(GET_SLOT(ctd, Matrix_xSym));              dtmp = REAL(GET_SLOT(ctd, Matrix_xSym));
# Line 370  Line 410
410              nnz = Lp[ncj];              nnz = Lp[ncj];
411              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, nnz));              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, nnz));
412              SET_SLOT(Ljj, Matrix_xSym, allocVector(REALSXP, nnz));              SET_SLOT(Ljj, Matrix_xSym, allocVector(REALSXP, nnz));
413              info = ldl_numeric(ncj, cp, ci, dtmp, Lp, parent, Lnz,              i = R_ldl_numeric(ncj, cp, ci, dtmp, Lp, parent,
414                                 INTEGER(GET_SLOT(Ljj, Matrix_iSym)),                                 INTEGER(GET_SLOT(Ljj, Matrix_iSym)),
415                                 REAL(GET_SLOT(Ljj, Matrix_xSym)),                                 REAL(GET_SLOT(Ljj, Matrix_xSym)),
416                                 (double *) R_alloc(ncj, sizeof(double)), /* D */                                 (double *) R_alloc(ncj, sizeof(double)), /* D */
(double *) R_alloc(ncj, sizeof(double)), /* Y */
(int *) R_alloc(ncj, sizeof(int)),       /* Pattern */
(int *) R_alloc(ncj, sizeof(int)),       /* Flag */
417                                 (int *) NULL, (int *) NULL);                                 (int *) NULL, (int *) NULL);
418              if (info < ncj) error("identity matrix is not positive definite?");              if (i < ncj) error("identity matrix is not positive definite?");
419              Free(Lnz); Free(tmp);              Free(iPerm);
420          } else {          } else {
421              for (i = 0; i <= ncj; i++) {              for (i = 0; i <= ncj; i++) {
422                  Lp[i] = 0;                  Lp[i] = 0;
# Line 396  Line 433
433          SET_SLOT(Linvj, Matrix_xSym, allocVector(REALSXP, nnz));          SET_SLOT(Linvj, Matrix_xSym, allocVector(REALSXP, nnz));
434          dtmp = REAL(GET_SLOT(Linvj, Matrix_xSym));          dtmp = REAL(GET_SLOT(Linvj, Matrix_xSym));
435          for (i = 0; i < nnz; i++) dtmp[i] = 1.;          for (i = 0; i < nnz; i++) dtmp[i] = 1.;
436  /* FIXME: Update any blocks below the diagonal in this column if necessary*/          for (k = j+1; k < nf; k++) { /* Update other blocks in this column */
437          for (k = j+1; k < nf; k++) { /* update remaining columns, if any */              symbolic_right_unit_mm(Linvj, VECTOR_ELT(ctab, Lind(k,j)));
438            }
439            for (k = j+1; k < nf; k++) { /* Update remaining columns */
440              for (i = k; i < nf; i++) block_update(ctab, j, k, i);              for (i = k; i < nf; i++) block_update(ctab, j, k, i);
441          }          }
442          for (i= 0; i < j; i++) { /* copy blocks to the left */          for (i= 0; i < j; i++) { /* copy blocks to the left */

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