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