SCM Repository

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

Diff of /pkg/src/bCrosstab.c

revision 414, Mon Jan 3 18:55:11 2005 UTC revision 415, Mon Jan 3 19:05:11 2005 UTC
# Line 120  Line 120
120      return -1;      return -1;
121  }  }
122
123    static int*
124    expand_column_pointers(int ncol, const int mp[], int mj[])
125    {
126        int j;
127        for (j = 0; j < ncol; j++) {
128            int j2 = mp[j+1], jj;
129            for (jj = mp[j]; jj < j2; jj++) mj[jj] = j;
130        }
131        return mj;
132    }
133
134  /**  /**
135   * Replace the structure of C by the structure of CA   * Replace the structure of C by the structure of CA
136   *   *
# Line 144  Line 155
155          cnr, cnz = length(cip),          cnr, cnz = length(cip),
156          i, j;          i, j;
157
158      if ((length(cpp) - 1) != anc)      if ((length(cpp) - 1) != anc) /* A is square so can compare no of cols */
159          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)",
160                anc, length(cpp) - 1);                anc, length(cpp) - 1);
161      if (anz < 1) return;        /* C is the identity */      if (anz < 1) return;        /* A is the identity */
162      cnr = -1;                   /* number of rows in C is max(ci + 1) */      cnr = -1;                   /* number of rows in C is max(ci + 1) */
163      for (i = 0; i < cnz; i++) {      for (i = 0; i < cnz; i++) {
164          int ri = ci[i] + 1;          int ri = ci[i] + 1;
# Line 201  Line 212
212  }  }
213
214  /**  /**
215     * Replace the structure of C by the structure of CA'
216     *
217     * @param A a unit lower triangular cscBlocked object
218     * @param C a cscBlocked object to be updated
219     */
220    static void
221    symbolic_right_unit_mm_trans(SEXP A, SEXP C)
222    {
223        SEXP aip = GET_SLOT(A, Matrix_iSym),
224            app = GET_SLOT(A, Matrix_pSym),
225            cip = GET_SLOT(C, Matrix_iSym),
226            cpp = GET_SLOT(C, Matrix_pSym);
227        int *ai = INTEGER(aip),
228            *ap = INTEGER(app),
229            *ci = INTEGER(cip),
230            *cp = INTEGER(cpp),
231            anc = length(app) - 1,
232            anz = length(aip),
233            cnz = length(cip),
234            j, nextra = 0;
235
236        if ((length(cpp) - 1) != anc)
237            error("No. of cols in A (%d) does not match no. of cols in C (%d)",
238                  anc, length(cpp) - 1);
239        if (anz < 1) return;        /* A is the identity */
240        for (j = 0; j < anc; j++) { /* bound the number of extra triplets */
241            int aj2 = ap[j + 1], cj2 = cp[j + 1], ka, kc;
242            for (ka = ap[j]; ka < aj2; ka++) {
243                for (kc = cp[j]; kc < cj2; kc++) {
244                    if (check_csc_index(cp, ci, ai[ka], ci[kc]) < 0) nextra++;
245                }
246            }
247        }
248        if (nextra) {
249            int cnr, ntot = cnz + nextra, pos;
250            int *dims = INTEGER(getAttrib(GET_SLOT(C, Matrix_xSym), R_DimSymbol)),
251                *Ti = Memcpy((int *) Calloc(ntot, int), ci, cnz),
252                *Tj = expand_column_pointers(anc, cp, Calloc(ntot, int)),
253                *Ci = Calloc(ntot, int);
254
255            for (j = 0, pos = cnz; j < anc; j++) {
256                int aj2 = ap[j + 1], cj2 = cp[j + 1], ka, kc;
257                for (ka = ap[j]; ka < aj2; ka++) {
258                    for (kc = cp[j]; kc < cj2; kc++) {
259                        if (check_csc_index(cp, ci, ci[kc], ai[ka]) < 0) {
260                            Tj[pos] = ai[ka];
261                            Ti[pos] = ci[kc];
262                            pos++;
263                        }
264                    }
265                }
266            }
267            for (j = 0, cnr = -1; j < cnz; j++) {int rr = ci[j]; if (rr > cnr) cnr = rr;}
268            cnr++;                  /* maximum index is 1 less than number of rows */
269            triplet_to_col(cnr, anc, ntot, Ti, Tj, (double *) NULL,
270                           INTEGER(cpp), Ci, (double *) NULL);
271            cnz = cp[anc];
272            SET_SLOT(C, Matrix_iSym, allocVector(INTSXP, cnz));
273            SET_SLOT(C, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], cnz));
274            Free(Ti); Free(Tj); Free(Ci);
275        }
276    }
277
278    /**
279   * Update a block of L in the blocked crosstabulation   * Update a block of L in the blocked crosstabulation
280   *   *
281   * @param ctab pointer to a blocked crosstabulation object   * @param ctab pointer to a blocked crosstabulation object
# Line 289  Line 364
364      }      }
365  }  }
366
static int*
expand_column_pointers(int ncol, int nnz, 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;
}

367  /**  /**
368   * Permute the levels of one of the grouping factors in a bCrosstab object   * Permute the levels of one of the grouping factors in a bCrosstab object
369   *   *
# Line 324  Line 388
388          int *cp = INTEGER(GET_SLOT(cscb, Matrix_pSym)),          int *cp = INTEGER(GET_SLOT(cscb, Matrix_pSym)),
389              nnz = length(cscbi);              nnz = length(cscbi);
390          double *cx = REAL(GET_SLOT(cscb, Matrix_xSym));          double *cx = REAL(GET_SLOT(cscb, Matrix_xSym));
391          int *mj = expand_column_pointers(ncol, nnz, cp, Calloc(nnz, int));          int *mj = expand_column_pointers(ncol, cp, Calloc(nnz, int));
392          int *mi = Memcpy(Calloc(nnz, int), INTEGER(cscbi), nnz);          int *mi = Memcpy(Calloc(nnz, int), INTEGER(cscbi), nnz);
393          double *mx = Memcpy(Calloc(nnz, double), cx, nnz);          double *mx = Memcpy(Calloc(nnz, double), cx, nnz);
394
# Line 497  Line 561
561          dtmp = REAL(GET_SLOT(Linvj, Matrix_xSym));          dtmp = REAL(GET_SLOT(Linvj, Matrix_xSym));
562          for (i = 0; i < nnz; i++) dtmp[i] = 1.;          for (i = 0; i < nnz; i++) dtmp[i] = 1.;
563          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 */
564              symbolic_right_unit_mm(Linvj, VECTOR_ELT(ZZpO, Lind(k,j)));              symbolic_right_unit_mm_trans(Linvj, VECTOR_ELT(ZZpO, Lind(k,j)));
565          }          }
566          for (k = j+1; k < nf; k++) { /* Update remaining columns */          for (k = j+1; k < nf; k++) { /* Update remaining columns */
567              for (i = k; i < nf; i++) block_update(ZZpO, j, k, i);              for (i = k; i < nf; i++) block_update(ZZpO, j, k, i);

Legend:
 Removed from v.414 changed lines Added in v.415