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

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