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 447 - (view) (download) (as text)

1 : bates 347 #include "bCrosstab.h"
2 : bates 410 /* TODO:
3 : bates 431 * - Use the off-diagonal blocks of L
4 :     * - Remove the off-diagonal blocks of ZZpO
5 :     * - Only do a fill-reducing permutation on the first non-nested factor
6 : bates 410 */
7 : bates 347
8 :     /**
9 : bates 410 * Allocate a 3-dimensional array
10 :     *
11 :     * @param TYP The R Type code (e.g. INTSXP)
12 :     * @param nr number of rows
13 :     * @param nc number of columns
14 :     * @param nf number of faces
15 :     *
16 :     * @return A 3-dimensional array of the indicated dimensions and type
17 :     */
18 :     static
19 :     SEXP alloc3Darray(int TYP, int nr, int nc, int nf)
20 :     {
21 :     SEXP val, dd = PROTECT(allocVector(INTSXP, 3));
22 :    
23 :     INTEGER(dd)[0] = nr; INTEGER(dd)[1] = nc; INTEGER(dd)[2] = nf;
24 :     val = allocArray(TYP, dd);
25 :     UNPROTECT(1);
26 :     return val;
27 :     }
28 :    
29 :     /**
30 : bates 347 * Calculate the zero-based index in a row-wise packed lower
31 :     * triangular matrix. This is used for the arrays of blocked sparse matrices.
32 :     *
33 :     * @param i row number (0-based)
34 :     * @param k column number (0-based)
35 :     *
36 :     * @return The 0-based index of the (i,k) element of a row-wise packed lower
37 :     * triangular matrix.
38 :     */
39 : bates 370 static R_INLINE int
40 :     Lind(int i, int k)
41 : bates 347 {
42 :     return (i * (i + 1))/2 + k;
43 :     }
44 :    
45 :     /**
46 : bates 370 * Apply a permutation to an index vector
47 : bates 347 *
48 :     * @param i vector of 0-based indices
49 :     * @param nnz length of vector i
50 :     * @param perm 0-based permutation vector of length max(i)
51 :     */
52 : bates 370 static R_INLINE void
53 :     ind_permute(int i[], int nnz, const int perm[])
54 : bates 347 {
55 :     int j;
56 :     for (j = 0; j < nnz; j++) i[j] = perm[i[j]];
57 :     }
58 :    
59 :     /**
60 : bates 356 * Force indices to be in the upper triangle of a matrix
61 : bates 347 *
62 :     * @param i vector of 0-based row indices
63 :     * @param j vector of 0-based column indices
64 :     * @param nnz length of index vectors
65 :     */
66 : bates 370 static R_INLINE void
67 :     make_upper_triangular(int i[], int j[], int nnz)
68 : bates 347 {
69 :     int k;
70 :     for (k = 0; k < nnz; k++) {
71 : bates 356 if (i[k] > j[k]) {
72 : bates 347 int tmp = i[k];
73 :     i[k] = j[k];
74 :     j[k] = tmp;
75 :     }
76 :     }
77 :     }
78 :    
79 :     /**
80 : bates 410 * Create a named vector of type TYP
81 : bates 347 *
82 : bates 410 * @param TYP a vector SEXP type (e.g. REALSXP)
83 :     * @param names names of list elements with null string appended
84 : bates 347 *
85 : bates 410 * @return pointer to a named vector of type TYP
86 : bates 347 */
87 : bates 370 static SEXP
88 : bates 410 make_named(int TYP, char **names)
89 : bates 347 {
90 : bates 410 SEXP ans, nms;
91 :     int i, n;
92 : bates 347
93 : bates 410 for (n = 0; strlen(names[n]) > 0; n++) {}
94 :     ans = PROTECT(allocVector(TYP, n));
95 :     nms = PROTECT(allocVector(STRSXP, n));
96 : bates 356 for (i = 0; i < n; i++) SET_STRING_ELT(nms, i, mkChar(names[i]));
97 : bates 347 setAttrib(ans, R_NamesSymbol, nms);
98 :     UNPROTECT(2);
99 :     return ans;
100 :     }
101 :    
102 : bates 370 /**
103 :     * Check for the existence of the (row, col) pair in a csc structure.
104 :     *
105 :     * @param p vector of column pointers
106 :     * @param i vector of row indices
107 :     * @param row row index
108 :     * @param col column index
109 :     *
110 :     * @return index of element at (row, col) if it exists, otherwise -1
111 :     */
112 :     static R_INLINE int
113 :     check_csc_index(const int p[], const int i[], int row, int col)
114 : bates 356 {
115 :     int k, k2 = p[col + 1];
116 :     /* use a linear search for now */
117 :     /* perhaps replace by bsearch later */
118 :     for (k = p[col]; k < k2; k++) {
119 :     if (i[k] == row) return k;
120 :     }
121 :     return -1;
122 :     }
123 : bates 370
124 : bates 415 static int*
125 :     expand_column_pointers(int ncol, const int mp[], int mj[])
126 :     {
127 :     int j;
128 :     for (j = 0; j < ncol; j++) {
129 :     int j2 = mp[j+1], jj;
130 :     for (jj = mp[j]; jj < j2; jj++) mj[jj] = j;
131 :     }
132 :     return mj;
133 :     }
134 :    
135 : bates 347 /**
136 : bates 442 * 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 : bates 347 *
139 : bates 442 * @param anc number of columns in A
140 :     * @param Parent parent array for A
141 : bates 410 * @param C a cscBlocked object to be updated
142 : bates 347 */
143 : bates 370 static void
144 : bates 442 symbolic_right_unit_sm(int anc, const int Parent[], SEXP C)
145 : bates 347 {
146 : bates 442 SEXP cip = GET_SLOT(C, Matrix_iSym),
147 : bates 370 cpp = GET_SLOT(C, Matrix_pSym);
148 :     int *Flag,
149 :     *ci = INTEGER(cip),
150 :     *cp = INTEGER(cpp),
151 :     *ncp,
152 :     cnr, cnz = length(cip),
153 :     i, j;
154 :    
155 : bates 415 if ((length(cpp) - 1) != anc) /* A is square so can compare no of cols */
156 : bates 370 error("No. of rows in A (%d) does not match no. of cols in C (%d)",
157 : bates 415 anc, length(cpp) - 1);
158 : bates 442 i = 1; /* check for A being the identity */
159 :     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 : bates 370 for (i = 0; i < cnz; i++) {
168 :     int ri = ci[i] + 1;
169 :     if (cnr < ri) cnr = ri;
170 :     }
171 :     Flag = Calloc(cnr, int);
172 :     ncp = Calloc(anc + 1, int); /* new column pointers */
173 :    
174 :     ncp[0] = 0;
175 :     for (j = 0; j < anc; j++) {
176 : bates 442 int cj2 = cp[j + 1], kk, kc;
177 : bates 370 for (i = 0; i < cnr; i++) Flag[i] = 0;
178 :     ncp[j+1] = ncp[j] + cj2 - cp[j];
179 :     /* positions of current column j of C */
180 :     for (kc = cp[j]; kc < cj2; kc++) Flag[ci[kc]] = 1;
181 :     /* other positions in column j of product */
182 : bates 442 for (kk = Parent[j]; kk >= 0; kk = Parent[kk]) { /* loop over the Parent */
183 :     int kk2 = cp[kk + 1];
184 : bates 370 for (kc = cp[kk]; kc < kk2; kc++) {
185 :     if (!Flag[ci[kc]]) {
186 :     ncp[j+1]++;
187 :     Flag[ci[kc]] = 1;
188 :     }
189 : bates 356 }
190 :     }
191 :     }
192 : bates 370 if (ncp[anc] > cp[anc]) {
193 : bates 410 int *dims, *nci, nnz = ncp[anc], pos = 0;
194 : bates 370 double *ncx;
195 : bates 356
196 : bates 370 SET_SLOT(C, Matrix_iSym, allocVector(INTSXP, nnz));
197 :     nci = INTEGER(GET_SLOT(C, Matrix_iSym));
198 : bates 410 dims = INTEGER(getAttrib(GET_SLOT(C, Matrix_xSym), R_DimSymbol));
199 :     SET_SLOT(C, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], nnz));
200 : bates 370 ncx = REAL(GET_SLOT(C, Matrix_xSym));
201 :     for (i = 0; i < nnz; i++) ncx[i] = 1.;
202 :     /* As Diana Krall said, "Just do it again." */
203 :     for (j = 0; j < anc; j++) {
204 : bates 442 int cj2 = cp[j + 1], kc, kk;
205 : bates 370 for (i = 0; i < cnr; i++) Flag[i] = 0;
206 :     for (kc = cp[j]; kc < cj2; kc++) Flag[ci[kc]] = 1;
207 : bates 442 for (kk = Parent[j]; kk >= 0; kk = Parent[kk]) {
208 :     int kk2 = cp[kk + 1];
209 : bates 370 for (kc = cp[kk]; kc < kk2; kc++) Flag[ci[kc]] = 1;
210 : bates 356 }
211 : bates 370 for (i = 0; i < cnr; i++) if (Flag[i]) nci[pos++] = i;
212 : bates 356 }
213 : bates 370 Memcpy(cp, ncp, anc + 1);
214 :     }
215 :     Free(Flag); Free(ncp);
216 : bates 347 }
217 : bates 415
218 :     /**
219 : bates 442 * Replace the structure of C by the structure of CA^{-T}
220 : bates 415 *
221 :     * @param A a unit lower triangular cscBlocked object
222 :     * @param C a cscBlocked object to be updated
223 :     */
224 :     static void
225 : bates 442 symbolic_right_unit_mm_trans(int anc, const int Parent[], SEXP C)
226 : bates 415 {
227 : bates 442 SEXP cip = GET_SLOT(C, Matrix_iSym),
228 : bates 415 cpp = GET_SLOT(C, Matrix_pSym);
229 : bates 442 int *ci = INTEGER(cip),
230 : bates 415 *cp = INTEGER(cpp),
231 :     cnz = length(cip),
232 : bates 442 i, j, nextra = 0;
233 : bates 415
234 :     if ((length(cpp) - 1) != anc)
235 :     error("No. of cols in A (%d) does not match no. of cols in C (%d)",
236 :     anc, length(cpp) - 1);
237 : bates 442
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 : bates 415 for (j = 0; j < anc; j++) { /* bound the number of extra triplets */
248 : bates 442 int cj2 = cp[j + 1], ka, kc;
249 :     for (ka = Parent[j]; ka >= 0; ka = Parent[ka]) {
250 : bates 415 for (kc = cp[j]; kc < cj2; kc++) {
251 : bates 442 if (check_csc_index(cp, ci, ci[kc], ka) < 0) nextra++;
252 : bates 415 }
253 :     }
254 :     }
255 :     if (nextra) {
256 :     int cnr, ntot = cnz + nextra, pos;
257 :     int *dims = INTEGER(getAttrib(GET_SLOT(C, Matrix_xSym), R_DimSymbol)),
258 :     *Ti = Memcpy((int *) Calloc(ntot, int), ci, cnz),
259 :     *Tj = expand_column_pointers(anc, cp, Calloc(ntot, int)),
260 :     *Ci = Calloc(ntot, int);
261 :    
262 :     for (j = 0, pos = cnz; j < anc; j++) {
263 : bates 442 int cj2 = cp[j + 1], ka, kc;
264 :     for (ka = Parent[j]; ka >= 0; ka = Parent[ka]) {
265 : bates 415 for (kc = cp[j]; kc < cj2; kc++) {
266 : bates 442 if (check_csc_index(cp, ci, ci[kc], ka) < 0) {
267 :     Tj[pos] = ka;
268 : bates 415 Ti[pos] = ci[kc];
269 :     pos++;
270 :     }
271 :     }
272 :     }
273 :     }
274 : bates 442 for (j = 0, cnr = 0; j < cnz; j++) { /* determine number of rows in C */
275 :     int rr = ci[j] + 1;
276 :     if (rr > cnr) cnr = rr;
277 :     }
278 : bates 415 triplet_to_col(cnr, anc, ntot, Ti, Tj, (double *) NULL,
279 :     INTEGER(cpp), Ci, (double *) NULL);
280 :     cnz = cp[anc];
281 :     SET_SLOT(C, Matrix_iSym, allocVector(INTSXP, cnz));
282 :     SET_SLOT(C, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], cnz));
283 :     Free(Ti); Free(Tj); Free(Ci);
284 :     }
285 :     }
286 : bates 370
287 : bates 347 /**
288 : bates 370 * Update a block of L in the blocked crosstabulation
289 : bates 347 *
290 :     * @param ctab pointer to a blocked crosstabulation object
291 : bates 356 * @param j index of updating column
292 :     * @param k column index of block to be updated
293 :     * @param i row index of block to be updated (j < k <= i)
294 : bates 347 */
295 : bates 370 static void
296 : bates 447 block_update(SEXP L, SEXP ZZpO, int j, int k, int i)
297 : bates 347 {
298 : bates 447 SEXP tb = (i == k) ? VECTOR_ELT(ZZpO, i) : VECTOR_ELT(L, Lind(i, k)),
299 :     ib = VECTOR_ELT(L, Lind(i, j)),
300 :     kb = VECTOR_ELT(L, Lind(k, j));
301 : bates 356 SEXP tpp = GET_SLOT(tb, Matrix_pSym),
302 :     kpp = GET_SLOT(kb, Matrix_pSym);
303 :     int *ti = INTEGER(GET_SLOT(tb, Matrix_iSym)),
304 :     *tp = INTEGER(tpp),
305 :     *ii = INTEGER(GET_SLOT(ib, Matrix_iSym)),
306 :     *ip = INTEGER(GET_SLOT(ib, Matrix_pSym)),
307 :     *ki = INTEGER(GET_SLOT(kb, Matrix_iSym)),
308 :     *kp = INTEGER(kpp),
309 :     tnc = length(tpp) - 1,
310 :     knc = length(kpp) - 1;
311 :     int jj, extra;
312 : bates 347
313 : bates 356 if (k > i || j >= k)
314 :     error("i,j,k values of %d,%d,%d do not satisfy j < k <= i",
315 :     i, j, k);
316 :     /* bound the number of extra elements */
317 :     extra = 0;
318 :     for (jj = 0; jj < knc; jj++) {
319 :     int i1, kk, i2 = ip[jj + 1], k2 = kp[jj + 1];
320 :     for (kk = kp[jj]; kk < k2; kk++) {
321 :     for (i1 = ip[jj]; i1 < i2; i1++) {
322 :     if ((check_csc_index(tp, ti, ii[i1], ki[kk]) < 0) &&
323 :     /* only update upper triangle of
324 :     * diagonal blocks */
325 :     ((k != i) || (ii[i1] <= ki[kk]))) extra++;
326 :     }
327 :     }
328 : bates 347 }
329 : bates 356 if (!extra) return;
330 :     {
331 :     int pos, nnz = tp[tnc];
332 : bates 410 int ntot = nnz + extra, tnr;
333 : bates 356 int *Ai = Calloc(ntot, int),
334 :     *Ti = Calloc(ntot, int),
335 :     *Tj = Calloc(ntot, int),
336 : bates 410 *dims;
337 : bates 356 double *Ax;
338 :    
339 :     Memcpy(Ti, ti, nnz); /* make a copy of the row indices */
340 :     for (pos = 0, jj = 0; jj < tnc; jj++) { /* fill in the column indices */
341 :     int j2 = tp[jj + 1];
342 :     for (; pos < j2; pos++) Tj[pos] = jj;
343 :     }
344 :     /* add the extra elements */
345 :     for (jj = 0; jj < knc; jj++) {
346 :     int i1, kk, i2 = ip[jj + 1], k2 = kp[jj + 1];
347 :     for (kk = kp[jj]; kk < k2; kk++) {
348 :     for (i1 = ip[jj]; i1 < i2; i1++) {
349 :     if ((check_csc_index(tp, ti, ii[i1], ki[kk]) < 0) &&
350 :     ((k != i) || (ii[i1] <= ki[kk]))) {
351 :     Ti[pos] = ii[i1];
352 :     Tj[pos] = ki[kk];
353 :     pos++;
354 :     }
355 :     }
356 :     }
357 :     }
358 : bates 410 /* Pass nlev instead of doing this. The dimensions are nlev[i], nlev[k] */
359 :     /* Determine maximum row index in T */
360 :     tnr = -1; for (jj = 0; jj < ntot; jj++) if (Ti[jj] > tnr) tnr = Ti[jj];
361 :     tnr++; /* increment by 1 to get number of rows */
362 :     triplet_to_col(tnr, tnc, ntot, Ti, Tj, (double *) NULL,
363 : bates 356 tp, Ai, (double *) NULL);
364 :     nnz = tp[tnc];
365 :     SET_SLOT(tb, Matrix_iSym, allocVector(INTSXP, nnz));
366 :     Memcpy(INTEGER(GET_SLOT(tb, Matrix_iSym)), Ai, nnz);
367 : bates 410 dims = INTEGER(getAttrib(GET_SLOT(tb, Matrix_xSym), R_DimSymbol));
368 :     SET_SLOT(tb, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], nnz));
369 : bates 356 Ax = REAL(GET_SLOT(tb, Matrix_xSym));
370 :     for (j = 0; j < nnz; j++) Ax[j] = 1.;
371 :     Free(Ai); Free(Ti); Free(Tj);
372 :     return;
373 :     }
374 : bates 347 }
375 :    
376 :     /**
377 :     * Permute the levels of one of the grouping factors in a bCrosstab object
378 :     *
379 :     * @param ctab Pointer to a bCrosstab object
380 : bates 356 * @param nf number of factors in ctab
381 :     * @param jj index (0-based) of the factor levels to permute
382 :     * @param ncj number of columns in level jj
383 : bates 347 * @param perm permutation (0-based) to apply
384 : bates 356 * @param pperm inverse of the permutation
385 : bates 347 */
386 : bates 370 static void
387 : bates 447 bCrosstab_permute(SEXP ctab, int nf, int jj, const int nlev[], const int iperm[])
388 : bates 347 {
389 : bates 410 int j;
390 : bates 356 for (j = 0; j < nf; j++) {
391 : bates 410 int ind = (j < jj ? Lind(jj, j) : Lind(j, jj)),
392 :     ncol = (j < jj ? nlev[j] : nlev[jj]),
393 :     nrow = (j < jj ? nlev[jj] : nlev[j]);
394 :     SEXP cscb = VECTOR_ELT(ctab, ind),
395 :     cscbi = GET_SLOT(cscb, Matrix_iSym);
396 :     int *cp = INTEGER(GET_SLOT(cscb, Matrix_pSym)),
397 :     nnz = length(cscbi);
398 :     double *cx = REAL(GET_SLOT(cscb, Matrix_xSym));
399 : bates 415 int *mj = expand_column_pointers(ncol, cp, Calloc(nnz, int));
400 : bates 410 int *mi = Memcpy(Calloc(nnz, int), INTEGER(cscbi), nnz);
401 :     double *mx = Memcpy(Calloc(nnz, double), cx, nnz);
402 : bates 356
403 : bates 410 if (j <= jj) ind_permute(mi, nnz, iperm);
404 :     if (j >= jj) ind_permute(mj, nnz, iperm);
405 :     if (j == jj) make_upper_triangular(mi, mj, nnz);
406 :     triplet_to_col(nrow, ncol, nnz, mi, mj, mx, cp, INTEGER(cscbi), cx);
407 :     Free(mi); Free(mj); Free(mx);
408 : bates 347 }
409 :     }
410 :    
411 : bates 447 static void
412 :     symmetric_permute(SEXP A, int nlev, const int iperm[])
413 :     {
414 :     SEXP AiP = GET_SLOT(A, Matrix_iSym);
415 :     int *Ap = INTEGER(GET_SLOT(A, Matrix_pSym)),
416 :     nnz = length(AiP);
417 :     double *Ax = REAL(GET_SLOT(A, Matrix_xSym));
418 :     int *mj = expand_column_pointers(nlev, Ap, Calloc(nnz, int));
419 :     int *mi = Memcpy(Calloc(nnz, int), INTEGER(AiP), nnz);
420 :     double *mx = Memcpy(Calloc(nnz, double), Ax, nnz);
421 :    
422 :     ind_permute(mi, nnz, iperm);
423 :     ind_permute(mj, nnz, iperm);
424 :     make_upper_triangular(mi, mj, nnz);
425 :     triplet_to_col(nlev, nlev, nnz, mi, mj, mx, Ap, INTEGER(AiP), Ax);
426 :     Free(mi); Free(mj); Free(mx);
427 :     }
428 :    
429 : bates 370 /**
430 :     * Apply a permutation vector to the levels of a factor.
431 :     *
432 :     * The dest pointer is assumed to point to a copy of the src pointer's
433 :     * contents.
434 :     *
435 :     * @param dest pointer to the destination factor
436 :     * @param src pointer to the source factor
437 :     * @param perm permutation vector (0-based)
438 :     * @param iperm inverse permutation vector (0-based)
439 :     */
440 : bates 356 static void
441 :     factor_levels_permute(SEXP dest, SEXP src, const int perm[],
442 :     const int iperm[])
443 :     {
444 :     SEXP dlev = getAttrib(dest, R_LevelsSymbol),
445 :     slev = getAttrib(src, R_LevelsSymbol);
446 :     int nlev = length(dlev), flen = length(dest);
447 :     int *d = INTEGER(dest), *s = INTEGER(src), i;
448 :    
449 :     if (length(slev) != nlev)
450 :     error("number of levels in src and dest must match");
451 :     if (length(src) != flen)
452 :     error("length of src and dest must match");
453 :     for (i = 0; i < nlev; i++)
454 :     SET_STRING_ELT(dlev, i, STRING_ELT(slev, perm[i]));
455 :     for (i = 0; i < flen; i++)
456 :     d[i] = 1 + iperm[s[i]-1];
457 :     }
458 :    
459 : bates 370 /**
460 : bates 410 * Create and populate slots in an lmer object from the blocked crosstabulation.
461 : bates 370 *
462 : bates 410 * @param val Pointer to an lmer object
463 : bates 370 */
464 : bates 410 void
465 :     lmer_populate(SEXP val)
466 : bates 356 {
467 : bates 444 SEXP D, L, Parent, ZZpO, flist = GET_SLOT(val, Matrix_flistSym),
468 :     perm, Omega, ZtZ = GET_SLOT(val, Matrix_ZtZSym);
469 :     int j, k, nf = length(flist);
470 : bates 410 int *nc = INTEGER(GET_SLOT(val, Matrix_ncSym)), *Gp,
471 :     *nlev = Calloc(nf, int), npairs = (nf * (nf + 1))/2;
472 :     char *statnms[] = {"factored", "inverted", ""},
473 : bates 442 *devnms[] = {"ML", "REML", ""},
474 :     *pnms[] = {"index", "block", ""};
475 : bates 410
476 :     /* Allocate fixed-sized slots */
477 :     SET_SLOT(val, Matrix_statusSym, make_named(LGLSXP, statnms ));
478 :     SET_SLOT(val, Matrix_devianceSym, make_named(REALSXP, devnms));
479 :     SET_SLOT(val, Matrix_devCompSym, allocVector(REALSXP, 4));
480 : bates 444 /* Allocate slots that are lists of length nf */
481 :     SET_SLOT(val, Matrix_ZZpOSym, allocVector(VECSXP, nf));
482 : bates 410 ZZpO = GET_SLOT(val, Matrix_ZZpOSym);
483 : bates 444 setAttrib(ZZpO, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
484 : bates 410 SET_SLOT(val, Matrix_DSym, allocVector(VECSXP, nf));
485 :     D = GET_SLOT(val, Matrix_DSym);
486 :     setAttrib(D, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
487 :     SET_SLOT(val, Matrix_permSym, allocVector(VECSXP, nf));
488 :     perm = GET_SLOT(val, Matrix_permSym);
489 : bates 356 setAttrib(perm, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
490 : bates 410 SET_SLOT(val, Matrix_ParentSym, allocVector(VECSXP, nf));
491 :     Parent = GET_SLOT(val, Matrix_ParentSym);
492 : bates 356 setAttrib(Parent, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
493 : bates 410 SET_SLOT(val, Matrix_OmegaSym, allocVector(VECSXP, nf));
494 :     Omega = GET_SLOT(val, Matrix_OmegaSym);
495 :     setAttrib(Omega, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
496 :    
497 :     /* Allocate peculiar length slots */
498 :     SET_SLOT(val, Matrix_LSym, allocVector(VECSXP, npairs));
499 :     L = GET_SLOT(val, Matrix_LSym);
500 :     SET_SLOT(val, Matrix_GpSym, allocVector(INTSXP, nf + 1));
501 :     Gp = INTEGER(GET_SLOT(val, Matrix_GpSym));
502 :     Gp[0] = 0;
503 : bates 356 for (j = 0; j < nf; j++) {
504 : bates 410 nlev[j] = length(getAttrib(VECTOR_ELT(flist, j), R_LevelsSymbol));
505 :     Gp[j + 1] = Gp[j] + nc[j] * nlev[j];
506 :     SET_VECTOR_ELT(D, j, alloc3Darray(REALSXP, nc[j], nc[j], nlev[j]));
507 :     SET_VECTOR_ELT(Omega, j, allocMatrix(REALSXP, nc[j], nc[j]));
508 : bates 444 SET_VECTOR_ELT(ZZpO, j, duplicate(VECTOR_ELT(ZtZ, Lind(j, j))));
509 :     for (k = j; k < nf; k++)
510 :     SET_VECTOR_ELT(L, Lind(k, j), duplicate(VECTOR_ELT(ZtZ, Lind(k, j))));
511 : bates 410 }
512 :     SET_SLOT(val, Matrix_XtXSym, allocMatrix(REALSXP, nc[nf], nc[nf]));
513 : bates 431 AZERO(REAL(GET_SLOT(val, Matrix_XtXSym)), nc[nf] * nc[nf]);
514 : bates 410 SET_SLOT(val, Matrix_RXXSym, allocMatrix(REALSXP, nc[nf], nc[nf]));
515 : bates 431 AZERO(REAL(GET_SLOT(val, Matrix_RXXSym)), nc[nf] * nc[nf]);
516 : bates 410 SET_SLOT(val, Matrix_ZtXSym, allocMatrix(REALSXP, Gp[nf], nc[nf]));
517 :     SET_SLOT(val, Matrix_RZXSym, allocMatrix(REALSXP, Gp[nf], nc[nf]));
518 :     for (j = 0; j < nf; j++) {
519 : bates 444 int dind = Lind(j, j), i;
520 :     SEXP ctd = VECTOR_ELT(ZZpO, j); /* diagonal in crosstab */
521 :     SEXP Ljj = VECTOR_ELT(L, dind),
522 :     cpp = GET_SLOT(ctd, Matrix_pSym),
523 : bates 442 cip = GET_SLOT(ctd, Matrix_iSym), parent;
524 : bates 444 int *Lp = INTEGER(GET_SLOT(Ljj, Matrix_pSym)), *Perm,
525 :     *cp = INTEGER(cpp),
526 : bates 442 *ci = INTEGER(cip), *dims,
527 : bates 356 ncj = length(cpp) - 1,
528 :     nnz = length(cip);
529 :    
530 : bates 442 SET_VECTOR_ELT(Parent, j, make_named(VECSXP, pnms));
531 :     parent = VECTOR_ELT(Parent, j);
532 :     SET_VECTOR_ELT(parent, 0, allocVector(INTSXP, ncj));
533 :     SET_VECTOR_ELT(parent, 1, allocVector(INTSXP, ncj));
534 : bates 356 SET_VECTOR_ELT(perm, j, allocVector(INTSXP, ncj));
535 :     Perm = INTEGER(VECTOR_ELT(perm, j));
536 : bates 410 dims = INTEGER(getAttrib(GET_SLOT(ctd, Matrix_xSym), R_DimSymbol));
537 : bates 356 if (nnz > ncj) { /* calculate fill-reducing permutation */
538 : bates 410 SEXP fac = VECTOR_ELT(flist, j);
539 :     SEXP fcp = PROTECT(duplicate(fac));
540 : bates 370 int *iPerm = Calloc(ncj, int);
541 : bates 356
542 : bates 370 ssc_metis_order(ncj, cp, ci, Perm, iPerm);
543 : bates 447 /* apply to the crosstabulation, L, and ZZpO */
544 :     bCrosstab_permute(ZtZ, nf, j, nlev, iPerm);
545 :     bCrosstab_permute(L, nf, j, nlev, iPerm);
546 :     symmetric_permute(VECTOR_ELT(ZZpO, j), nlev[j], iPerm);
547 : bates 356 /* apply to the factor */
548 : bates 410 factor_levels_permute(fac, fcp, Perm, iPerm);
549 : bates 356 /* symbolic analysis to get Parent */
550 : bates 442 R_ldl_symbolic(ncj, cp, ci, Lp, INTEGER(VECTOR_ELT(parent, 0)),
551 : bates 356 (int *) NULL, (int *) NULL);
552 : bates 442 for (i = 0; i < ncj; i++)
553 :     INTEGER(VECTOR_ELT(parent, 1))[i] =
554 :     (INTEGER(VECTOR_ELT(parent, 0))[i] < 0) ? -1 : j;
555 : bates 356 nnz = Lp[ncj];
556 :     SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, nnz));
557 : bates 410 SET_SLOT(Ljj, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], nnz));
558 :     Free(iPerm); UNPROTECT(1);
559 : bates 356 } else {
560 : bates 442 for (i = 0; i < ncj; i++) {
561 : bates 356 Lp[i] = 0;
562 : bates 442 INTEGER(VECTOR_ELT(parent,0))[i] = -1;
563 :     INTEGER(VECTOR_ELT(parent,1))[i] = -1;
564 : bates 356 Perm[i] = i;
565 :     }
566 : bates 442 Lp[ncj] = 0;
567 : bates 356 SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, 0));
568 : bates 410 SET_SLOT(Ljj, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], 0));
569 : bates 356 }
570 : bates 370 for (k = j+1; k < nf; k++) { /* Update other blocks in this column */
571 : bates 442 symbolic_right_unit_mm_trans(ncj, INTEGER(VECTOR_ELT(parent, 0)),
572 : bates 444 VECTOR_ELT(L, Lind(k,j)));
573 : bates 370 }
574 :     for (k = j+1; k < nf; k++) { /* Update remaining columns */
575 : bates 447 for (i = k; i < nf; i++) block_update(L, ZZpO, j, k, i);
576 : bates 356 }
577 :     }
578 : bates 442 /* Convert blockwise Parent arrays to extended Parent arrays */
579 :     for (j = 0; j < (nf - 1); j++) { /* Parent[nf] does not need conversion */
580 :     SEXP Ljp1j = VECTOR_ELT(L, Lind(j + 1, j)),
581 :     LpP = GET_SLOT(Ljp1j, Matrix_pSym);
582 :     int *Li = INTEGER(GET_SLOT(Ljp1j, Matrix_iSym)),
583 :     *Lp = INTEGER(LpP),
584 :     *block = INTEGER(VECTOR_ELT(VECTOR_ELT(Parent, j), 1)),
585 :     *parent = INTEGER(VECTOR_ELT(VECTOR_ELT(Parent, j), 0)),
586 :     i, nlev = length(LpP) - 1;
587 :     for (i = 0; i < nlev; i++) {
588 :     if (block[i] < 0) {
589 :     block[i] = j + 1;
590 :     parent[i] = Li[Lp[i]];
591 :     }
592 :     }
593 :     }
594 : bates 410 Free(nlev);
595 : bates 356 }

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