SCM

SCM Repository

[matrix] Annotation of /pkg/src/bCrosstab.c
ViewVC logotype

Annotation of /pkg/src/bCrosstab.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 370 - (view) (download) (as text)

1 : bates 347 #include "bCrosstab.h"
2 :    
3 :     /**
4 :     * Calculate the zero-based index in a row-wise packed lower
5 :     * triangular matrix. This is used for the arrays of blocked sparse matrices.
6 :     *
7 :     * @param i row number (0-based)
8 :     * @param k column number (0-based)
9 :     *
10 :     * @return The 0-based index of the (i,k) element of a row-wise packed lower
11 :     * triangular matrix.
12 :     */
13 : bates 370 static R_INLINE int
14 :     Lind(int i, int k)
15 : bates 347 {
16 :     return (i * (i + 1))/2 + k;
17 :     }
18 :    
19 :     /**
20 : bates 370 * Apply a permutation to an index vector
21 : bates 347 *
22 :     * @param i vector of 0-based indices
23 :     * @param nnz length of vector i
24 :     * @param perm 0-based permutation vector of length max(i)
25 :     */
26 : bates 370 static R_INLINE void
27 :     ind_permute(int i[], int nnz, const int perm[])
28 : bates 347 {
29 :     int j;
30 :     for (j = 0; j < nnz; j++) i[j] = perm[i[j]];
31 :     }
32 :    
33 :     /**
34 : bates 356 * Force indices to be in the upper triangle of a matrix
35 : bates 347 *
36 :     * @param i vector of 0-based row indices
37 :     * @param j vector of 0-based column indices
38 :     * @param nnz length of index vectors
39 :     */
40 : bates 370 static R_INLINE void
41 :     make_upper_triangular(int i[], int j[], int nnz)
42 : bates 347 {
43 :     int k;
44 :     for (k = 0; k < nnz; k++) {
45 : bates 356 if (i[k] > j[k]) {
46 : bates 347 int tmp = i[k];
47 :     i[k] = j[k];
48 :     j[k] = tmp;
49 :     }
50 :     }
51 :     }
52 :    
53 :     /**
54 :     * Create a named list of length n
55 :     *
56 :     * @param n length of list to return
57 :     * @param names names of list elements
58 :     *
59 :     * @return pointer to a named list of length n
60 :     */
61 : bates 370 static SEXP
62 :     make_named_list(int n, char **names)
63 : bates 347 {
64 :     SEXP ans = PROTECT(allocVector(VECSXP, n)),
65 : bates 356 nms = PROTECT(allocVector(STRSXP, n));
66 : bates 347 int i;
67 :    
68 : bates 356 for (i = 0; i < n; i++) SET_STRING_ELT(nms, i, mkChar(names[i]));
69 : bates 347 setAttrib(ans, R_NamesSymbol, nms);
70 :     UNPROTECT(2);
71 :     return ans;
72 :     }
73 :    
74 : bates 370 /**
75 :     * 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 : bates 356 {
87 :     int k, k2 = p[col + 1];
88 :     /* use a linear search for now */
89 :     /* perhaps replace by bsearch later */
90 :     for (k = p[col]; k < k2; k++) {
91 :     if (i[k] == row) return k;
92 :     }
93 :     return -1;
94 :     }
95 : bates 370
96 : bates 347 /**
97 : bates 370 * Replace the structure of C by the structure of CA
98 : bates 347 *
99 : bates 370 * @param A a unit lower tscMatrix
100 :     * @param C cscMatrix to be updated
101 : bates 347 */
102 : bates 370 static void
103 :     symbolic_right_unit_mm(SEXP A, SEXP C)
104 : bates 347 {
105 : bates 370 SEXP aip = GET_SLOT(A, Matrix_iSym),
106 :     app = GET_SLOT(A, Matrix_pSym),
107 :     cip = GET_SLOT(C, Matrix_iSym),
108 :     cpp = GET_SLOT(C, Matrix_pSym);
109 :     int *Flag,
110 :     *ai = INTEGER(aip),
111 :     *ap = INTEGER(app),
112 :     *ci = INTEGER(cip),
113 :     *cp = INTEGER(cpp),
114 :     *ncp,
115 :     anc = length(app) - 1,
116 :     anz = length(aip),
117 :     cnr, cnz = length(cip),
118 :     i, j;
119 :    
120 :     if ((length(cpp) - 1) != anc)
121 :     error("No. of rows in A (%d) does not match no. of cols in C (%d)",
122 :     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 : bates 356 }
148 :     }
149 :     }
150 : bates 370 if (ncp[anc] > cp[anc]) {
151 :     int *nci, nnz = ncp[anc], pos = 0;
152 :     double *ncx;
153 : bates 356
154 : bates 370 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 : bates 356 }
168 : bates 370 for (i = 0; i < cnr; i++) if (Flag[i]) nci[pos++] = i;
169 : bates 356 }
170 : bates 370 Memcpy(cp, ncp, anc + 1);
171 :     }
172 :     Free(Flag); Free(ncp);
173 : bates 347 }
174 : bates 370
175 : bates 347 /**
176 : bates 370 * Update a block of L in the blocked crosstabulation
177 : bates 347 *
178 :     * @param ctab pointer to a blocked crosstabulation object
179 : bates 356 * @param j index of updating column
180 :     * @param k column index of block to be updated
181 :     * @param i row index of block to be updated (j < k <= i)
182 : bates 347 */
183 : bates 370 static void
184 :     block_update(SEXP ctab, int j, int k, int i)
185 : bates 347 {
186 : bates 356 SEXP tb = VECTOR_ELT(ctab, Lind(i, k)),
187 :     ib = VECTOR_ELT(ctab, Lind(i, j)),
188 :     kb = VECTOR_ELT(ctab, Lind(k, j));
189 :     SEXP tpp = GET_SLOT(tb, Matrix_pSym),
190 :     kpp = GET_SLOT(kb, Matrix_pSym);
191 :     int *ti = INTEGER(GET_SLOT(tb, Matrix_iSym)),
192 :     *tp = INTEGER(tpp),
193 :     *ii = INTEGER(GET_SLOT(ib, Matrix_iSym)),
194 :     *ip = INTEGER(GET_SLOT(ib, Matrix_pSym)),
195 :     *ki = INTEGER(GET_SLOT(kb, Matrix_iSym)),
196 :     *kp = INTEGER(kpp),
197 :     tnc = length(tpp) - 1,
198 :     knc = length(kpp) - 1;
199 :     int jj, extra;
200 : bates 347
201 : bates 356 if (k > i || j >= k)
202 :     error("i,j,k values of %d,%d,%d do not satisfy j < k <= i",
203 :     i, j, k);
204 :     /* bound the number of extra elements */
205 :     extra = 0;
206 :     for (jj = 0; jj < knc; jj++) {
207 :     int i1, kk, i2 = ip[jj + 1], k2 = kp[jj + 1];
208 :     for (kk = kp[jj]; kk < k2; kk++) {
209 :     for (i1 = ip[jj]; i1 < i2; i1++) {
210 :     if ((check_csc_index(tp, ti, ii[i1], ki[kk]) < 0) &&
211 :     /* only update upper triangle of
212 :     * diagonal blocks */
213 :     ((k != i) || (ii[i1] <= ki[kk]))) extra++;
214 :     }
215 :     }
216 : bates 347 }
217 : bates 356 if (!extra) return;
218 :     {
219 :     int pos, nnz = tp[tnc];
220 :     int ntot = nnz + extra;
221 :     int *Ai = Calloc(ntot, int),
222 :     *Ti = Calloc(ntot, int),
223 :     *Tj = Calloc(ntot, int),
224 :     *Dims = INTEGER(GET_SLOT(tb, Matrix_DimSym));
225 :     double *Ax;
226 :    
227 :     Memcpy(Ti, ti, nnz); /* make a copy of the row indices */
228 :     for (pos = 0, jj = 0; jj < tnc; jj++) { /* fill in the column indices */
229 :     int j2 = tp[jj + 1];
230 :     for (; pos < j2; pos++) Tj[pos] = jj;
231 :     }
232 :     /* add the extra elements */
233 :     for (jj = 0; jj < knc; jj++) {
234 :     int i1, kk, i2 = ip[jj + 1], k2 = kp[jj + 1];
235 :     for (kk = kp[jj]; kk < k2; kk++) {
236 :     for (i1 = ip[jj]; i1 < i2; i1++) {
237 :     if ((check_csc_index(tp, ti, ii[i1], ki[kk]) < 0) &&
238 :     ((k != i) || (ii[i1] <= ki[kk]))) {
239 :     Ti[pos] = ii[i1];
240 :     Tj[pos] = ki[kk];
241 :     pos++;
242 :     }
243 :     }
244 :     }
245 :     }
246 :     triplet_to_col(Dims[0], Dims[1], ntot, Ti, Tj, (double *) NULL,
247 :     tp, Ai, (double *) NULL);
248 :     nnz = tp[tnc];
249 :     SET_SLOT(tb, Matrix_iSym, allocVector(INTSXP, nnz));
250 :     Memcpy(INTEGER(GET_SLOT(tb, Matrix_iSym)), Ai, nnz);
251 :     SET_SLOT(tb, Matrix_xSym, allocVector(REALSXP, nnz));
252 :     Ax = REAL(GET_SLOT(tb, Matrix_xSym));
253 :     for (j = 0; j < nnz; j++) Ax[j] = 1.;
254 :     Free(Ai); Free(Ti); Free(Tj);
255 :     return;
256 :     }
257 : bates 347 }
258 :    
259 :     /**
260 :     * Permute the levels of one of the grouping factors in a bCrosstab object
261 :     *
262 :     * @param ctab Pointer to a bCrosstab object
263 : bates 356 * @param nf number of factors in ctab
264 :     * @param jj index (0-based) of the factor levels to permute
265 :     * @param ncj number of columns in level jj
266 : bates 347 * @param perm permutation (0-based) to apply
267 : bates 356 * @param pperm inverse of the permutation
268 : bates 347 */
269 : bates 370 static void
270 :     bCrosstab_permute(SEXP ctab, int nf, int jj,
271 :     const int perm[], const int iperm[])
272 : bates 347 {
273 : bates 356 SEXP trip, ipt, jpt;
274 :     int j, nnz;
275 : bates 347
276 : bates 356 for (j = 0; j < nf; j++) {
277 :     int ind = (j < jj ? Lind(jj, j) : Lind(j, jj));
278 :     SEXP csc = VECTOR_ELT(ctab, ind);
279 :     int *Dims = INTEGER(GET_SLOT(csc, Matrix_DimSym));
280 :    
281 :     trip = csc_to_triplet(csc);
282 : bates 347 ipt = GET_SLOT(trip, Matrix_iSym);
283 : bates 356 nnz = length(ipt);
284 : bates 347 jpt = GET_SLOT(trip, Matrix_jSym);
285 : bates 356 if (j <= jj) ind_permute(INTEGER(ipt), nnz, iperm);
286 :     if (j >= jj) ind_permute(INTEGER(jpt), nnz, iperm);
287 :     if (j == jj)
288 :     make_upper_triangular(INTEGER(ipt), INTEGER(jpt), nnz);
289 :     triplet_to_col(Dims[0], Dims[1], nnz, INTEGER(ipt),
290 :     INTEGER(jpt),
291 :     REAL(GET_SLOT(trip, Matrix_xSym)),
292 :     INTEGER(GET_SLOT(csc, Matrix_pSym)),
293 :     INTEGER(GET_SLOT(csc, Matrix_iSym)),
294 :     REAL(GET_SLOT(csc, Matrix_xSym)));
295 : bates 347 }
296 :     }
297 :    
298 : bates 370 /**
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 : bates 356 static void
310 :     factor_levels_permute(SEXP dest, SEXP src, const int perm[],
311 :     const int iperm[])
312 :     {
313 :     SEXP dlev = getAttrib(dest, R_LevelsSymbol),
314 :     slev = getAttrib(src, R_LevelsSymbol);
315 :     int nlev = length(dlev), flen = length(dest);
316 :     int *d = INTEGER(dest), *s = INTEGER(src), i;
317 :    
318 :     if (length(slev) != nlev)
319 :     error("number of levels in src and dest must match");
320 :     if (length(src) != flen)
321 :     error("length of src and dest must match");
322 :     for (i = 0; i < nlev; i++)
323 :     SET_STRING_ELT(dlev, i, STRING_ELT(slev, perm[i]));
324 :     for (i = 0; i < flen; i++)
325 :     d[i] = 1 + iperm[s[i]-1];
326 :     }
327 :    
328 : bates 370 /**
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 : bates 356 {
342 :     char *anms[] = {"flist", "L", "Linv", "perm", "Parent"};
343 :     SEXP flist = VECTOR_ELT(bCtab, 0),
344 :     ctab = PROTECT(duplicate(VECTOR_ELT(bCtab, 1))),
345 : bates 370 ans = PROTECT(make_named_list(sizeof(anms)/sizeof(char *), anms));
346 : bates 356 SEXP fcp, perm, L, Linv, Parent;
347 :     int ctbl = length(ctab), j, nf = length(flist);
348 :    
349 :     if (ctbl != (nf*(nf + 1))/2)
350 :     error("length of ctab = %d is not permisable", ctbl);
351 :     SET_VECTOR_ELT(ans, 0, duplicate(flist));
352 :     SET_VECTOR_ELT(ans, 1, allocVector(VECSXP, ctbl));
353 :     SET_VECTOR_ELT(ans, 2, allocVector(VECSXP, nf));
354 :     SET_VECTOR_ELT(ans, 3, allocVector(VECSXP, nf));
355 :     SET_VECTOR_ELT(ans, 4, allocVector(VECSXP, nf));
356 :     fcp = VECTOR_ELT(ans, 0);
357 :     L = VECTOR_ELT(ans, 1);
358 :     setAttrib(L, R_NamesSymbol, duplicate(getAttrib(ctab, R_NamesSymbol)));
359 :     Linv = VECTOR_ELT(ans, 2);
360 :     setAttrib(Linv, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
361 :     perm = VECTOR_ELT(ans, 3);
362 :     setAttrib(perm, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
363 :     Parent = VECTOR_ELT(ans, 4);
364 :     setAttrib(Parent, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
365 :     for (j = 0; j < nf; j++) {
366 :     int dind = Lind(j, j), i, k;
367 :     SEXP ctd = VECTOR_ELT(ctab, dind); /* diagonal in crosstab */
368 :     SEXP Dimslot = GET_SLOT(ctd, Matrix_DimSym),
369 : bates 359 Linvj, Ljj, cpp = GET_SLOT(ctd, Matrix_pSym),
370 : bates 356 cip = GET_SLOT(ctd, Matrix_iSym);
371 : bates 359 int *Lp, *Linvp, *Perm, *cp = INTEGER(cpp),
372 : bates 356 *ci = INTEGER(cip), *parent,
373 :     ncj = length(cpp) - 1,
374 :     nnz = length(cip);
375 : bates 359 double *dtmp;
376 : bates 356
377 :     SET_VECTOR_ELT(Parent, j, allocVector(INTSXP, ncj));
378 :     parent = INTEGER(VECTOR_ELT(Parent, j));
379 :     SET_VECTOR_ELT(perm, j, allocVector(INTSXP, ncj));
380 :     Perm = INTEGER(VECTOR_ELT(perm, j));
381 :     SET_VECTOR_ELT(L, dind, NEW_OBJECT(MAKE_CLASS("cscMatrix")));
382 :     Ljj = VECTOR_ELT(L, dind);
383 : bates 359 SET_VECTOR_ELT(Linv, j, NEW_OBJECT(MAKE_CLASS("cscMatrix")));
384 :     Linvj = VECTOR_ELT(Linv, j);
385 : bates 356 SET_SLOT(Ljj, Matrix_DimSym, duplicate(Dimslot));
386 : bates 359 SET_SLOT(Linvj, Matrix_DimSym, duplicate(Dimslot));
387 : bates 356 SET_SLOT(Ljj, Matrix_factorization, allocVector(VECSXP, 0));
388 : bates 359 SET_SLOT(Linvj, Matrix_factorization, allocVector(VECSXP, 0));
389 : bates 356 SET_SLOT(Ljj, Matrix_pSym, allocVector(INTSXP, ncj + 1));
390 : bates 359 SET_SLOT(Linvj, Matrix_pSym, allocVector(INTSXP, ncj + 1));
391 : bates 356 Lp = INTEGER(GET_SLOT(Ljj, Matrix_pSym));
392 : bates 359 Linvp = INTEGER(GET_SLOT(Linvj, Matrix_pSym));
393 : bates 356 if (nnz > ncj) { /* calculate fill-reducing permutation */
394 : bates 370 int *iPerm = Calloc(ncj, int);
395 : bates 356
396 : bates 370 ssc_metis_order(ncj, cp, ci, Perm, iPerm);
397 : bates 356 /* apply to the crosstabulation */
398 : bates 370 bCrosstab_permute(ctab, nf, j, Perm, iPerm);
399 : bates 356 /* apply to the factor */
400 :     factor_levels_permute(VECTOR_ELT(fcp, j),
401 : bates 370 VECTOR_ELT(flist, j), Perm, iPerm);
402 : bates 356 /* symbolic analysis to get Parent */
403 : bates 370 R_ldl_symbolic(ncj, cp, ci, Lp, parent,
404 : bates 356 (int *) NULL, (int *) NULL);
405 : bates 359 /* decompose the identity to get the row pointers */
406 :     dtmp = REAL(GET_SLOT(ctd, Matrix_xSym));
407 :     for (i = 0; i < nnz; i++) dtmp[i] = 0.; /* initialize */
408 :     /* diagonal el is last element in the column */
409 :     for (i = 0; i < ncj; i++) dtmp[cp[i+1] - 1] = 1.;
410 : bates 356 nnz = Lp[ncj];
411 :     SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, nnz));
412 :     SET_SLOT(Ljj, Matrix_xSym, allocVector(REALSXP, nnz));
413 : bates 370 i = R_ldl_numeric(ncj, cp, ci, dtmp, Lp, parent,
414 : bates 359 INTEGER(GET_SLOT(Ljj, Matrix_iSym)),
415 :     REAL(GET_SLOT(Ljj, Matrix_xSym)),
416 :     (double *) R_alloc(ncj, sizeof(double)), /* D */
417 :     (int *) NULL, (int *) NULL);
418 : bates 370 if (i < ncj) error("identity matrix is not positive definite?");
419 :     Free(iPerm);
420 : bates 356 } else {
421 :     for (i = 0; i <= ncj; i++) {
422 :     Lp[i] = 0;
423 :     parent[i] = -1;
424 :     Perm[i] = i;
425 :     }
426 :     SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, 0));
427 :     SET_SLOT(Ljj, Matrix_xSym, allocVector(REALSXP, 0));
428 :     }
429 : bates 359 /* Derive the diagonal block of Linv */
430 :     nnz = parent_inv_ap(ncj, 0, parent, Linvp);
431 :     SET_SLOT(Linvj, Matrix_iSym, allocVector(INTSXP, nnz));
432 :     parent_inv_ai(ncj, 0, parent, INTEGER(GET_SLOT(Linvj, Matrix_iSym)));
433 :     SET_SLOT(Linvj, Matrix_xSym, allocVector(REALSXP, nnz));
434 :     dtmp = REAL(GET_SLOT(Linvj, Matrix_xSym));
435 :     for (i = 0; i < nnz; i++) dtmp[i] = 1.;
436 : bates 370 for (k = j+1; k < nf; k++) { /* Update other blocks in this column */
437 :     symbolic_right_unit_mm(Linvj, VECTOR_ELT(ctab, Lind(k,j)));
438 :     }
439 :     for (k = j+1; k < nf; k++) { /* Update remaining columns */
440 : bates 356 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 */
443 :     SET_VECTOR_ELT(L, Lind(j,i), duplicate(VECTOR_ELT(ctab, Lind(j,i))));
444 :     }
445 :     }
446 :     UNPROTECT(2);
447 :     return ans;
448 :     }
449 :    
450 :    
451 :    
452 : bates 347

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