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