SCM

SCM Repository

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

Annotation of /pkg/src/ssclme.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : bates 10 #include "ssclme.h"
2 :    
3 :     static
4 :     void ssclme_copy_ctab(int nf, const int nc[], SEXP ctab, SEXP ssc)
5 :     {
6 :     int *snc, i, copyonly = 1;
7 :    
8 :     for (i = 0; i < nf; i++) {
9 :     if (nc[i] > 1) copyonly = 0;
10 :     }
11 :     if (copyonly) {
12 :     SET_SLOT(ssc, Matrix_pSym, duplicate(GET_SLOT(ctab, Matrix_pSym)));
13 :     SET_SLOT(ssc, Matrix_iSym, duplicate(GET_SLOT(ctab, Matrix_iSym)));
14 :     SET_SLOT(ssc, Matrix_xSym, duplicate(GET_SLOT(ctab, Matrix_xSym)));
15 :     SET_SLOT(ssc, Matrix_DimSym,
16 :     duplicate(GET_SLOT(ctab, Matrix_DimSym)));
17 :     SET_SLOT(ssc, Matrix_GpSym, duplicate(GET_SLOT(ctab, Matrix_GpSym)));
18 :     } else {
19 :     int
20 :     *GpIn = INTEGER(GET_SLOT(ctab, Matrix_GpSym)),
21 :     *GpOut,
22 :     *AiIn = INTEGER(GET_SLOT(ctab, Matrix_iSym)),
23 :     *AiOut,
24 :     *ApIn = INTEGER(GET_SLOT(ctab, Matrix_pSym)),
25 :     *ApOut,
26 :     nIn = GpIn[nf], nOut, nzOut,
27 :     *dims,
28 :     *map = Calloc(nIn + 1, int), /* column number map */
29 :     *ncc = Calloc(nIn, int); /* number of columns out for this
30 :     * col in */
31 :    
32 :     SET_SLOT(ssc, Matrix_GpSym, allocVector(INTSXP, nf + 1));
33 :     GpOut = INTEGER(GET_SLOT(ssc, Matrix_GpSym));
34 :     map[0] = GpOut[0] = 0;
35 :     for (i = 0; i < nf; i++) {
36 :     int j, GpIni = GpIn[i], GpInip1 = GpIn[i+1], nci = nc[i];
37 :     GpOut[i+1] = GpOut[i] + (GpInip1 - GpIni) * nci;
38 :     for (j = GpIni; j < GpInip1; j++) {
39 :     ncc[j] = nci;
40 :     map[j+1] = map[j] + nci;
41 :     }
42 :     }
43 :     nOut = GpOut[nf]; /* size of output matrix */
44 :     SET_SLOT(ssc, Matrix_DimSym,
45 :     duplicate(GET_SLOT(ctab, Matrix_DimSym)));
46 :     dims = INTEGER(GET_SLOT(ssc, Matrix_DimSym));
47 :     dims[0] = dims[1] = nOut;
48 :     SET_SLOT(ssc, Matrix_pSym, allocVector(INTSXP, nOut + 1));
49 :     ApOut = INTEGER(GET_SLOT(ssc, Matrix_pSym));
50 :     ApOut[0] = 0;
51 :     for (i = 0; i < nf; i++) { /* determine the column pointers */
52 :     int j, jout = GpOut[i], nci = nc[i], p2 = GpIn[i+1];
53 :     for (j = GpIn[i]; j < p2; j++) {
54 :     int k, nj = 0, p3 = ApIn[j+1];
55 :     for (k = ApIn[j]; k < p3; k++) {
56 :     nj += ncc[AiIn[k]];
57 :     }
58 :     nj -= nci - 1;
59 :     ApOut[jout+1] = ApOut[jout] + nj;
60 :     jout++;
61 :     for (k = 1; k < nci; k++) {
62 :     ApOut[jout+1] = ApOut[jout] + nj + k;
63 :     jout++;
64 :     }
65 :     }
66 :     }
67 :     nzOut = ApOut[nOut]; /* number of non-zeros in output */
68 :     SET_SLOT(ssc, Matrix_xSym, allocVector(REALSXP, nzOut));
69 :     memset(REAL(GET_SLOT(ssc, Matrix_xSym)), 0,
70 :     sizeof(double) * nzOut);
71 :     SET_SLOT(ssc, Matrix_iSym, allocVector(INTSXP, nzOut));
72 :     AiOut = INTEGER(GET_SLOT(ssc, Matrix_iSym));
73 :     for (i = 0; i < nf; i++) { /* fill in the rows */
74 :     int j, jj, nci = nc[i], p2 = GpIn[i+1];
75 :     for (j = GpIn[i]; j < p2; j++) { /* first col for input col */
76 : bates 168 int k, mj = map[j], p3 = ApIn[j+1], pos = ApOut[mj];
77 :     for (k = ApIn[j]; k < p3; k++) {
78 :     int ii = AiIn[k], ncci = ncc[ii];
79 :     AiOut[pos++] = map[ii];
80 :     if (ii < j) { /* above the diagonal */
81 :     for (jj = 1; jj < ncci; jj++) {
82 :     AiOut[pos+1] = AiOut[pos] + 1;
83 :     pos++;
84 :     }
85 : bates 10 }
86 : bates 168 for (jj = 1; jj < nci; jj++) { /* repeat the column adding
87 :     * another diagonal element */
88 :     int mjj = mj + jj, pj = ApOut[mjj], pjm1 = ApOut[mjj-1];
89 :     Memcpy(AiOut + pj, AiOut + pjm1, pj - pjm1);
90 :     AiOut[ApOut[mjj + 1] - 1] = mjj; /* maybe mjj-1? */
91 :     }
92 : bates 10 }
93 :     }
94 :     }
95 :     Free(map); Free(ncc);
96 :     }
97 :     SET_SLOT(ssc, Matrix_ncSym, allocVector(INTSXP, nf + 2));
98 :     snc = INTEGER(GET_SLOT(ssc, Matrix_ncSym));
99 :     for (i = 0; i <= nf; i++) {
100 :     snc[i] = nc[i];
101 :     }
102 :     }
103 :    
104 : bates 159 /**
105 :     * Calculate and store the maximum number of off-diagonal elements in
106 :     * the inverse of L, based on the elimination tree. The maximum is
107 :     * itself stored in the Parent array. (FIXME: come up with a better design.)
108 :     *
109 :     * @param n number of columns in the matrix
110 :     * @param Parent elimination tree for the matrix
111 :     */
112 :     static void ssclme_calc_maxod(int n, int Parent[])
113 : bates 10 {
114 : bates 159 int *sz = Calloc(n, int), i, mm = -1;
115 : bates 10 for (i = n - 1; i >= 0; i--) {
116 : bates 159 sz[i] = (Parent[i] < 0) ? 0 : (1 + sz[Parent[i]]);
117 :     if (sz[i] > mm) mm = sz[i];
118 : bates 10 }
119 : bates 159 Parent[n] = mm;
120 : bates 10 Free(sz);
121 :     }
122 :    
123 :     SEXP
124 : bates 168 ssclme_create(SEXP facs, SEXP ncv)
125 : bates 10 {
126 :     SEXP ctab, nms, ssc, tmp,
127 : bates 21 val = PROTECT(allocVector(VECSXP, 2)),
128 :     dd = PROTECT(allocVector(INTSXP, 3)); /* dimensions of 3-D arrays */
129 : bates 159 int *Ai, *Ap, *Gp, *Lp, *Parent,
130 : bates 21 *nc, Lnz, i, nf = length(facs), nzcol, pp1,
131 :     *dims = INTEGER(dd);
132 : bates 10
133 :     if (length(ncv) != (nf + 1))
134 :     error("length of nc (%d) should be length of facs (%d) + 1",
135 :     length(ncv), nf);
136 :     SET_VECTOR_ELT(val, 0, NEW_OBJECT(MAKE_CLASS("ssclme")));
137 :     ssc = VECTOR_ELT(val, 0);
138 :     /* Pairwise cross-tabulation */
139 :     ctab = PROTECT(sscCrosstab(facs, ScalarLogical(1)));
140 : bates 143 SET_VECTOR_ELT(val, 1, sscCrosstab_groupedPerm(ctab));
141 : bates 10 if (length(VECTOR_ELT(val, 1)) > 0) {/* Fill-reducing permutation */
142 :     ssc_symbolic_permute(INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],
143 :     1, INTEGER(VECTOR_ELT(val, 1)),
144 :     INTEGER(GET_SLOT(ctab, Matrix_pSym)),
145 :     INTEGER(GET_SLOT(ctab, Matrix_iSym)));
146 :     }
147 :     ssclme_copy_ctab(nf, INTEGER(ncv), ctab, ssc);
148 :     UNPROTECT(1); /* ctab */
149 :    
150 :     nzcol = INTEGER(GET_SLOT(ssc, Matrix_DimSym))[1];
151 :     Gp = INTEGER(GET_SLOT(ssc, Matrix_GpSym));
152 :     Ap = INTEGER(GET_SLOT(ssc, Matrix_pSym));
153 :     Ai = INTEGER(GET_SLOT(ssc, Matrix_iSym));
154 :     nc = INTEGER(GET_SLOT(ssc, Matrix_ncSym));
155 :     nc[nf + 1] = length(VECTOR_ELT(facs, 0)); /* number of observations */
156 :     /* Create slots */
157 :     pp1 = nc[nf];
158 :     SET_SLOT(ssc, Matrix_XtXSym, allocMatrix(REALSXP, pp1, pp1));
159 :     SET_SLOT(ssc, Matrix_RXXSym, allocMatrix(REALSXP, pp1, pp1));
160 :     SET_SLOT(ssc, Matrix_ZtXSym, allocMatrix(REALSXP, nzcol, pp1));
161 :     SET_SLOT(ssc, Matrix_RZXSym, allocMatrix(REALSXP, nzcol, pp1));
162 : bates 101 /* Zero symmetric matrices (cosmetic) */
163 : bates 10 memset(REAL(GET_SLOT(ssc, Matrix_XtXSym)), 0,
164 :     sizeof(double) * pp1 * pp1);
165 :     memset(REAL(GET_SLOT(ssc, Matrix_RXXSym)), 0,
166 :     sizeof(double) * pp1 * pp1);
167 :     SET_SLOT(ssc, Matrix_LpSym, allocVector(INTSXP, nzcol + 1));
168 :     Lp = INTEGER(GET_SLOT(ssc, Matrix_LpSym));
169 : bates 159 SET_SLOT(ssc, Matrix_ParentSym, allocVector(INTSXP, nzcol + 1));
170 : bates 10 Parent = INTEGER(GET_SLOT(ssc, Matrix_ParentSym));
171 :     SET_SLOT(ssc, Matrix_DSym, allocVector(REALSXP, nzcol));
172 :     SET_SLOT(ssc, Matrix_DIsqrtSym, allocVector(REALSXP, nzcol));
173 :     ldl_symbolic(nzcol, Ap, Ai, Lp, Parent,
174 :     (int *) R_alloc(nzcol, sizeof(int)), /* Lnz */
175 :     (int *) R_alloc(nzcol, sizeof(int)), /* Flag */
176 :     (int *) NULL, (int *) NULL); /* P & Pinv */
177 : bates 159 ssclme_calc_maxod(nzcol, Parent);
178 : bates 10 Lnz = Lp[nzcol];
179 :     SET_SLOT(ssc, Matrix_LiSym, allocVector(INTSXP, Lnz));
180 :     SET_SLOT(ssc, Matrix_LxSym, allocVector(REALSXP, Lnz));
181 :     SET_SLOT(ssc, Matrix_OmegaSym, allocVector(VECSXP, nf));
182 :     tmp = GET_SLOT(ssc, Matrix_OmegaSym);
183 :     setAttrib(tmp, R_NamesSymbol, getAttrib(facs, R_NamesSymbol));
184 :     for (i = 0; i < nf; i++) {
185 :     SET_VECTOR_ELT(tmp, i, allocMatrix(REALSXP, nc[i], nc[i]));
186 :     memset(REAL(VECTOR_ELT(tmp, i)), 0,
187 :     sizeof(double) * nc[i] * nc[i]);
188 :     }
189 :     SET_SLOT(ssc, Matrix_devianceSym, allocVector(REALSXP, 2));
190 :     tmp = GET_SLOT(ssc, Matrix_devianceSym);
191 :     setAttrib(tmp, R_NamesSymbol, allocVector(STRSXP, 2));
192 :     nms = getAttrib(tmp, R_NamesSymbol);
193 :     SET_STRING_ELT(nms, 0, mkChar("ML"));
194 :     SET_STRING_ELT(nms, 1, mkChar("REML"));
195 :     SET_SLOT(ssc, Matrix_devCompSym, allocVector(REALSXP, 4));
196 :     SET_SLOT(ssc, Matrix_statusSym, allocVector(LGLSXP, 2));
197 :     tmp = GET_SLOT(ssc, Matrix_statusSym);
198 :     LOGICAL(tmp)[0] = LOGICAL(tmp)[1] = 0;
199 :     setAttrib(tmp, R_NamesSymbol, allocVector(STRSXP, 2));
200 :     nms = getAttrib(tmp, R_NamesSymbol);
201 :     SET_STRING_ELT(nms, 0, mkChar("factored"));
202 :     SET_STRING_ELT(nms, 1, mkChar("inverted"));
203 :     SET_SLOT(ssc, Matrix_bVarSym, allocVector(VECSXP, nf));
204 :     tmp = GET_SLOT(ssc, Matrix_bVarSym);
205 :     setAttrib(tmp, R_NamesSymbol, getAttrib(facs, R_NamesSymbol));
206 :     for (i = 0; i < nf; i++) {
207 : bates 21 int nci = nc[i], mi = (Gp[i+1] - Gp[i])/nc[i];
208 : bates 10
209 : bates 21 dims[0] = dims[1] = nci;
210 :     dims[2] = mi;
211 :     SET_VECTOR_ELT(tmp, i, allocArray(REALSXP, dd));
212 : bates 10 memset(REAL(VECTOR_ELT(tmp, i)), 0,
213 : bates 21 sizeof(double) * nci * nci * mi);
214 : bates 10 }
215 : bates 21 UNPROTECT(2);
216 : bates 10 return val;
217 :     }
218 :    
219 :     static
220 :     void bVj_to_A(int ncj, int Gpj, int Gpjp, const double bVj[],
221 :     const int Ap[], const int Ai[], double Ax[])
222 :     {
223 :     int i, diag, k;
224 :     for (i = Gpj; i < Gpjp; i += ncj) {
225 :     for (k = 0; k < ncj; k++) {
226 :     diag = Ap[i + k + 1] - 1;
227 :     if (Ai[diag] != i+k)
228 :     error("Expected Ai[%d] to be %d (i.e on diagonal) not %d",
229 :     diag, i+k, Ai[diag]);
230 :     Memcpy(Ax + diag - k, bVj + (i+k-Gpj)*ncj, k + 1);
231 :     }
232 :     }
233 :     }
234 :    
235 :     SEXP
236 : bates 22 ssclme_transfer_dimnames(SEXP x, SEXP facs, SEXP mmats)
237 :     {
238 :     SEXP bVar = GET_SLOT(x, Matrix_bVarSym),
239 :     nms2 = PROTECT(allocVector(VECSXP, 2)),
240 :     nms3 = PROTECT(allocVector(VECSXP, 3));
241 :     int i, nf = length(mmats) - 1;
242 :     SEXP xcols = VECTOR_ELT(GetArrayDimnames(VECTOR_ELT(mmats, nf)), 1);
243 :    
244 :     for (i = 0; i < nf; i++) {
245 :     SEXP cnms = VECTOR_ELT(GetArrayDimnames(VECTOR_ELT(mmats, i)), 1);
246 :     SET_VECTOR_ELT(nms3, 0, cnms);
247 :     SET_VECTOR_ELT(nms3, 1, cnms);
248 :     SET_VECTOR_ELT(nms3, 2,
249 :     getAttrib(VECTOR_ELT(facs, i), R_LevelsSymbol));
250 :     dimnamesgets(VECTOR_ELT(bVar, i), duplicate(nms3));
251 :     }
252 :     SET_VECTOR_ELT(nms2, 0, xcols);
253 :     SET_VECTOR_ELT(nms2, 1, xcols);
254 :     dimnamesgets(GET_SLOT(x, Matrix_XtXSym), nms2);
255 :     dimnamesgets(GET_SLOT(x, Matrix_RXXSym), nms2);
256 :     UNPROTECT(2);
257 :     return R_NilValue;
258 :     }
259 :    
260 :     SEXP
261 : bates 10 ssclme_update_mm(SEXP x, SEXP facs, SEXP mmats)
262 :     {
263 : bates 22 SEXP bVar = GET_SLOT(x, Matrix_bVarSym);
264 : bates 10 int
265 :     *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
266 :     *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
267 :     *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
268 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
269 : bates 143 *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
270 : bates 10 i, j, k,
271 :     ione = 1,
272 :     nf = length(mmats) - 1,
273 :     nobs = nc[nf + 1],
274 :     nzcol = Gp[nf],
275 :     nz = Ap[nzcol],
276 :     pp1 = nc[nf];
277 :     double
278 :     **Z = Calloc(nf + 1, double *),
279 :     *Ax = REAL(GET_SLOT(x, Matrix_xSym)),
280 :     *XtX = REAL(GET_SLOT(x, Matrix_XtXSym)),
281 :     *ZtX = REAL(GET_SLOT(x, Matrix_ZtXSym)),
282 :     one = 1.0,
283 :     zero = 0.0;
284 :    
285 :     for (i = 0; i <= nf; i++) {
286 :     int *dims = INTEGER(getAttrib(VECTOR_ELT(mmats, i), R_DimSymbol)),
287 :     nci = nc[i];
288 :     if (nobs != dims[0])
289 :     error("Expected %d rows in the %d'th model matrix. Got %d",
290 :     nobs, i+1, dims[0]);
291 :     if (nci != dims[1])
292 :     error("Expected %d columns in the %d'th model matrix. Got %d",
293 :     nci, i+1, dims[1]);
294 :     Z[i] = REAL(VECTOR_ELT(mmats, i));
295 :     }
296 :     /* Create XtX - X is Z[nf] */
297 :     F77_CALL(dsyrk)("U", "T", nc+nf, &nobs, &one,
298 :     Z[nf], &nobs, &zero, XtX, nc + nf);
299 :     /* Zero the accumulators */
300 :     memset((void *) ZtX, 0, sizeof(double) * nzcol * pp1);
301 :     memset((void *) Ax, 0, sizeof(double) * nz);
302 :     for (j = 0; j < nf; j++) { /* Create ZtX */
303 :     int *fpj = INTEGER(VECTOR_ELT(facs, j)), ncj = nc[j],
304 :     Ncj = ncj > 1;
305 :     double
306 :     *bVj = REAL(VECTOR_ELT(bVar, j)),
307 :     *Zj = Z[j],
308 :     *zxj = ZtX + Gp[j];
309 :    
310 :     if (Ncj) { /* bVj will accumulate Z'Z blocks */
311 :     memset(bVj, 0, sizeof(double) * ncj * (Gp[j+1]-Gp[j]));
312 :     }
313 :     for (i = 0; i < nobs; i++) { /* accumulate diagonal of ZtZ */
314 :     int fpji = fpj[i] - 1, /* factor indices are 1-based */
315 :     dind = Ap[Gp[j] + fpji * ncj + 1] - 1;
316 :     if (Ai[dind] != (Gp[j] + fpji * ncj))
317 :     error("logic error in ssclme_update_mm");
318 :     if (Ncj) { /* use bVar to accumulate */
319 :     F77_CALL(dsyrk)("U", "T", &ncj, &ione, &one, Zj+i,
320 :     &nobs, &one, bVj + fpji*ncj*ncj, &ncj);
321 :     } else { /* update scalars directly */
322 :     Ax[dind] += Zj[i] * Zj[i];
323 :     }
324 :     /* update rows of Z'X */
325 :     F77_CALL(dgemm)("T", "N", &ncj, &pp1, &ione, &one,
326 :     Zj + i, &nobs, Z[nf] + i, &nobs,
327 :     &one, zxj + fpji * ncj, &nzcol);
328 :     }
329 :     if (Ncj) bVj_to_A(ncj, Gp[j], Gp[j+1], bVj, Ap, Ai, Ax);
330 :     for (k = j+1; k < nf; k++) { /* off-diagonals */
331 :     int *fpk = INTEGER(VECTOR_ELT(facs, k)),
332 :     *Apk = Ap + Gp[k],
333 : bates 111 nck = nc[k],
334 :     scalar = ncj == 1 && nck == 1;
335 : bates 10 double
336 : bates 111 *Zk = Z[k], *work;
337 :     if (!scalar) work = Calloc(ncj * nck, double);
338 : bates 10 for (i = 0; i < nobs; i++) {
339 :     int ii, ind = -1, fpji = fpj[i] - 1,
340 :     row = Gp[j] + fpji * ncj,
341 :     fpki = fpk[i] - 1,
342 : bates 111 lastind = Apk[fpki*nck + 1];
343 :     for (ii = Apk[fpki*nck]; ii < lastind; ii++) {
344 : bates 10 if (Ai[ii] == row) {
345 :     ind = ii;
346 :     break;
347 :     }
348 :     }
349 :     if (ind < 0) error("logic error in ssclme_update_mm");
350 : bates 111 if (scalar) { /* update scalars directly */
351 : bates 108 Ax[ind] += Zj[i] * Zk[i];
352 : bates 111 } else {
353 :     int jj, offset = ind - Apk[fpki * nck];
354 :     F77_CALL(dgemm)("T", "N", &ncj, &nck, &ione, &one,
355 :     Zj + i, &nobs, Zk + i, &nobs,
356 :     &zero, work, &ncj);
357 :     for (jj = 0; jj < nck; jj++) {
358 :     ind = Apk[fpki * nck + jj] + offset;
359 :     if (Ai[ind] != row)
360 :     error("logic error in ssclme_update_mm");
361 :     for (ii = 0; ii < ncj; ii++) {
362 :     Ax[ind++] += work[jj * ncj + ii];
363 :     }
364 :     }
365 : bates 10 }
366 :     }
367 : bates 111 if (!scalar) Free(work);
368 : bates 10 }
369 :     }
370 :     Free(Z);
371 : bates 22 ssclme_transfer_dimnames(x, facs, mmats);
372 : bates 143 status[0] = status[1] = 0;
373 : bates 10 return R_NilValue;
374 :     }
375 :    
376 : bates 97 SEXP ssclme_inflate_and_factor(SEXP x)
377 : bates 10 {
378 :     SEXP
379 : bates 97 GpSlot = GET_SLOT(x, Matrix_GpSym),
380 :     Omega = GET_SLOT(x, Matrix_OmegaSym);
381 :     int n = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];
382 : bates 10 int
383 : bates 97 *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
384 :     *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
385 : bates 10 *Flag = Calloc(n, int),
386 :     *Gp = INTEGER(GpSlot),
387 :     *Lnz = Calloc(n, int),
388 :     *Pattern = Calloc(n, int),
389 : bates 97 *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
390 : bates 10 j,
391 :     nf = length(GpSlot) - 1;
392 :     double
393 : bates 97 *D = REAL(GET_SLOT(x, Matrix_DSym)),
394 :     *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
395 : bates 10 *Y = Calloc(n, double),
396 :     *xcp = Calloc(Ap[n], double);
397 :    
398 : bates 97 Memcpy(xcp, REAL(GET_SLOT(x, Matrix_xSym)), Ap[n]);
399 : bates 10 for (j = 0; j < nf; j++) {
400 :     int diag, i, ii, k, G2 = Gp[j + 1], ncj = nc[j];
401 :     double *omgj = REAL(VECTOR_ELT(Omega, j));
402 :    
403 :     for (i = Gp[j]; i < G2; i += ncj) {
404 :     for (k = 0; k < ncj; k++) {
405 :     diag = Ap[i + k + 1] - 1;
406 :     if (Ai[diag] != i+k)
407 :     error("Expected Ai[%d] to be %d (i.e on diagonal) not %d",
408 :     diag, i+k, Ai[diag]);
409 :     for (ii = 0; ii <= k; ii++) {
410 :     xcp[diag + ii - k] += omgj[k*ncj + ii];
411 :     }
412 :     }
413 :     }
414 :     }
415 :     j = ldl_numeric(n, Ap, Ai, xcp,
416 : bates 97 INTEGER(GET_SLOT(x, Matrix_LpSym)),
417 :     INTEGER(GET_SLOT(x, Matrix_ParentSym)),
418 :     Lnz, INTEGER(GET_SLOT(x, Matrix_LiSym)),
419 :     REAL(GET_SLOT(x, Matrix_LxSym)),
420 : bates 10 D, Y, Pattern, Flag,
421 :     (int *) NULL, (int *) NULL); /* P & Pinv */
422 :     if (j != n)
423 :     error("rank deficiency of ZtZ+W detected at column %d",
424 :     j + 1);
425 :     for (j = 0; j < n; j++) DIsqrt[j] = 1./sqrt(D[j]);
426 :     Free(Lnz); Free(Flag); Free(Pattern); Free(Y); Free(xcp);
427 :     return R_NilValue;
428 :     }
429 :    
430 : bates 97 SEXP ssclme_factor(SEXP x)
431 : bates 10 {
432 : bates 97 int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
433 : bates 10
434 :     if (!status[0]) {
435 :     SEXP
436 : bates 97 GpSlot = GET_SLOT(x, Matrix_GpSym),
437 :     Omega = GET_SLOT(x, Matrix_OmegaSym);
438 : bates 10 int
439 :     *Gp = INTEGER(GpSlot),
440 : bates 97 *Li = INTEGER(GET_SLOT(x, Matrix_LiSym)),
441 :     *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
442 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
443 : bates 10 i,
444 : bates 97 n = INTEGER(GET_SLOT(x, Matrix_DimSym))[1],
445 : bates 10 nf = length(GpSlot) - 1,
446 :     nobs = nc[nf + 1],
447 :     nreml = nobs + 1 - nc[nf],
448 :     pp1 = nc[nf],
449 :     pp2 = pp1 + 1;
450 :     double
451 : bates 97 *D = REAL(GET_SLOT(x, Matrix_DSym)),
452 :     *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
453 :     *Lx = REAL(GET_SLOT(x, Matrix_LxSym)),
454 :     *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
455 :     *RZX = REAL(GET_SLOT(x, Matrix_RZXSym)),
456 :     *dcmp = REAL(getAttrib(x, Matrix_devCompSym)),
457 :     *deviance = REAL(getAttrib(x, Matrix_devianceSym)),
458 : bates 10 minus1 = -1.,
459 :     one = 1.;
460 :    
461 : bates 97 ssclme_inflate_and_factor(x);
462 : bates 10 /* Accumulate logdet of ZtZ+W */
463 :     dcmp[0] = dcmp[1] = dcmp[2] = dcmp[3] = 0.;
464 :     for (i = 0; i < n; i++) dcmp[0] += log(D[i]);
465 :     /* Accumulate logdet of W */
466 :     for (i = 0; i < nf; i++) {
467 :     int nci = nc[i],
468 :     mi = (Gp[i+1] - Gp[i])/nci;
469 :    
470 :     if (nci < 2) {
471 :     dcmp[1] += mi * log(REAL(VECTOR_ELT(Omega, i))[0]);
472 :     } else {
473 :     int j;
474 :     double
475 :     *tmp = Calloc(nci * nci, double),
476 :     accum = 0.;
477 :     F77_CALL(dpotrf)("U", &nci,
478 :     Memcpy(tmp, REAL(VECTOR_ELT(Omega, i)),
479 :     nci * nci),
480 :     &nci, &j);
481 :     if (j)
482 :     error("Omega[%d] is not positive definite", i + 1);
483 :     for (j = 0; j < nci; j++) {
484 :     accum += 2 * log(tmp[j * (nci + 1)]);
485 :     }
486 :     dcmp[1] += mi * accum;
487 :     Free(tmp);
488 :     }
489 :     }
490 :     /* ldl_lsolve on Z'X */
491 : bates 97 Memcpy(RZX, REAL(GET_SLOT(x, Matrix_ZtXSym)), n * pp1);
492 : bates 10 for (i = 0; i < pp1; i++) {
493 :     int j;
494 :     double *RZXi = RZX + i * n;
495 :     ldl_lsolve(n, RZXi, Lp, Li, Lx);
496 :     for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
497 :     }
498 :     /* downdate and factor X'X */
499 : bates 97 Memcpy(RXX, REAL(GET_SLOT(x, Matrix_XtXSym)), pp1 * pp1);
500 : bates 10 F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,
501 :     RZX, &n, &one, RXX, &pp1);
502 :     F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);
503 : bates 122 if (i) {
504 :     warning("Could not factor downdated X'X, code %d", i);
505 :     dcmp[2] = dcmp[3] = deviance[0] = deviance[1] = NA_REAL;
506 :     } else {
507 : bates 10 /* logdet of RXX */
508 : bates 122 for (i = 0; i < (pp1 - 1); i++)
509 :     dcmp[2] += 2 * log(RXX[i*pp2]);
510 : bates 10 /* logdet of Ryy */
511 : bates 122 dcmp[3] = 2. * log(RXX[pp1 * pp1 - 1]);
512 :     deviance[0] = /* ML criterion */
513 :     dcmp[0] - dcmp[1] + nobs*(1+dcmp[3]+log(2*PI/nobs));
514 :     deviance[1] = dcmp[0] - dcmp[1] + /* REML */
515 :     dcmp[2] + nreml * (1. + dcmp[3] + log(2. * PI/nreml));
516 :     }
517 : bates 10 status[0] = 1;
518 :     status[1] = 0;
519 :     }
520 :     return R_NilValue;
521 :     }
522 :    
523 :     static
524 :     int ldl_update_ind(int probe, int start, const int ind[])
525 :     {
526 :     while (ind[start] < probe) start++;
527 :     if (ind[start] > probe) error("logic error in ldl_inverse");
528 :     return start;
529 :     }
530 :    
531 :     /**
532 : bates 159 * Update the diagonal blocks of the inverse of LDL' (=Z'Z+W)
533 : bates 10 *
534 :     * @param x pointer to an ssclme object
535 :     *
536 :     * @return R_NilValue (x is updated in place)
537 :    
538 :     */
539 : bates 97 static
540 : bates 10 SEXP ldl_inverse(SEXP x)
541 :     {
542 :     SEXP
543 :     Gpsl = GET_SLOT(x, Matrix_GpSym),
544 :     bVar = GET_SLOT(x, Matrix_bVarSym);
545 :     int *Gp = INTEGER(Gpsl),
546 : bates 159 *Parent = INTEGER(GET_SLOT(x, Matrix_ParentSym)),
547 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
548 : bates 10 i,
549 :     nf = length(Gpsl) - 1,
550 : bates 159 nzc = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];
551 :     int maxod = Parent[nzc];
552 :     double *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym));
553 : bates 10
554 :     ssclme_factor(x);
555 : bates 159 if (maxod == 0) { /* L and L^{-1} are the identity */
556 : bates 10 for (i = 0; i < nf; i++) {
557 :     Memcpy(REAL(VECTOR_ELT(bVar, i)), DIsqrt + Gp[i],
558 :     Gp[i+1] - Gp[i]);
559 :     }
560 : bates 159 } else {
561 :     int *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
562 :     *Li = INTEGER(GET_SLOT(x, Matrix_LiSym));
563 :    
564 :     double one = 1.0, zero = 0.,
565 :     *Lx = REAL(GET_SLOT(x, Matrix_LxSym));
566 : bates 157
567 : bates 10 for (i = 0; i < nf; i++) {
568 : bates 158 int j, jj, k, kk, nci = nc[i], nr, p, p2, pj, pp,
569 : bates 157 m = maxod + 1,
570 : bates 159 *ind = Calloc(m, int), G1 = Gp[i], G2 = Gp[i+1];
571 : bates 10 double
572 :     *tmp = Calloc(m * nci, double),
573 : bates 159 *bVi = REAL(VECTOR_ELT(bVar, i));
574 : bates 10
575 : bates 159 /* initialize bVi to zero */
576 :     memset(bVi, 0, sizeof(double) * (G2 - G1) * nci);
577 :    
578 :     for (j = G1; j < G2; j += nci) {
579 : bates 158 kk = 0; /* ind gets indices of non-zeros */
580 : bates 10 jj = j; /* in this block of columns */
581 :     while (jj >= 0) {
582 :     ind[kk++] = jj;
583 :     jj = Parent[jj];
584 :     }
585 :     nr = kk; /* number of non-zeros in this block */
586 :     while (kk < m) ind[kk++] = nzc; /* placeholders */
587 :    
588 :     for (k = 0; k < nci; k++) {
589 :     double *ccol = tmp + k * nr;
590 :    
591 : bates 157 for (kk = 0; kk < nr; kk++) ccol[kk] = 0.;
592 :     ccol[k] = 1.; /* initialize from unit diagonal */
593 : bates 10 for (jj = j + k; jj >= 0; jj = Parent[jj]) {
594 :     p2 = Lp[jj+1];
595 : bates 158 pp = pj = ldl_update_ind(jj, 0, ind);
596 : bates 10 for (p = Lp[jj]; p < p2; p++) {
597 :     pp = ldl_update_ind(Li[p], pp, ind);
598 : bates 158 ccol[pp] -= Lx[p] * ccol[pj];
599 : bates 10 }
600 :     }
601 :     }
602 : bates 157
603 :     for (kk = 0; kk < nr; kk++) { /* scale rows */
604 : bates 10 for (k = 0; k < nci; k++) {
605 :     tmp[k * nr + kk] *= DIsqrt[ind[kk]];
606 :     }
607 :     }
608 : bates 113 F77_CALL(dsyrk)("L", "T", &nci, &nr, &one, tmp, &nr,
609 : bates 159 &zero, bVi + (j - G1)*nci, &nci);
610 :     F77_CALL(dpotrf)("L", &nci, bVi + (j - G1)*nci,
611 :     &nci, &jj);
612 :     if (jj) /* should never happen */
613 : bates 10 error(
614 : bates 159 "Rank deficient variance matrix at group %d, level %d, error code %d",
615 :     i + 1, j + 1, jj);
616 : bates 10 }
617 :     Free(tmp); Free(ind);
618 :     }
619 :     }
620 :     return R_NilValue;
621 :     }
622 : bates 101
623 : bates 97 SEXP ssclme_invert(SEXP x)
624 : bates 10 {
625 : bates 97 int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
626 :     if (!status[0]) ssclme_factor(x);
627 : bates 122 if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0]))
628 :     error("Unable to invert singular factor of downdated X'X");
629 : bates 10 if (!status[1]) {
630 :     SEXP
631 : bates 97 RZXsl = GET_SLOT(x, Matrix_RZXSym);
632 : bates 10 int
633 :     *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
634 : bates 97 *Li = INTEGER(GET_SLOT(x, Matrix_LiSym)),
635 :     *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
636 : bates 10 i,
637 :     n = dims[0],
638 :     pp1 = dims[1];
639 :     double
640 : bates 97 *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
641 :     *Lx = REAL(GET_SLOT(x, Matrix_LxSym)),
642 :     *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
643 : bates 10 *RZX = REAL(RZXsl),
644 :     one = 1.;
645 :    
646 :     F77_CALL(dtrtri)("U", "N", &pp1, RXX, &pp1, &i);
647 :     if (i)
648 :     error("DTRTRI returned error code %d", i);
649 :     F77_CALL(dtrmm)("R", "U", "N", "N", &n, &pp1, &one,
650 :     RXX, &pp1, RZX, &n);
651 :     for (i = 0; i < pp1; i++) {
652 :     int j; double *RZXi = RZX + i * n;
653 :     for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
654 :     ldl_ltsolve(n, RZXi, Lp, Li, Lx);
655 :     }
656 : bates 97 ldl_inverse(x);
657 : bates 10 status[1] = 1;
658 :     }
659 :     return R_NilValue;
660 :     }
661 :    
662 :     SEXP ssclme_initial(SEXP x)
663 :     {
664 :     SEXP Gpsl = GET_SLOT(x, Matrix_GpSym),
665 :     Omg = GET_SLOT(x, Matrix_OmegaSym);
666 :     int *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
667 :     *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
668 :     *Gp = INTEGER(Gpsl),
669 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
670 :     *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
671 :     i, nf = length(Gpsl) - 1;
672 :    
673 :     double *Ax = REAL(GET_SLOT(x, Matrix_xSym));
674 :    
675 :     for (i = 0; i < nf; i++) {
676 :     int
677 :     Gpi = Gp[i],
678 :     j, k,
679 :     nci = nc[i],
680 :     ncip1 = nci + 1,
681 :     p2 = Gp[i+1];
682 :     double
683 :     mi = 0.375 * ((double) nci)/((double) (p2 - Gpi)),
684 :     *mm = REAL(VECTOR_ELT(Omg, i));
685 :    
686 :     memset((void *) mm, 0, sizeof(double) * nci * nci);
687 :     for (j = Gpi; j < p2; j += nci) {
688 :     for (k = 0; k < nci; k++) {
689 :     int jk = j+k, jj = Ap[jk+1] - 1;
690 :     if (Ai[jj] != jk) error("malformed ZtZ structure");
691 :     mm[k * ncip1] += Ax[jj] * mi;
692 :     }
693 :     }
694 :     }
695 :     status[0] = status[1] = 0;
696 :     return R_NilValue;
697 :     }
698 :    
699 :     /**
700 :     * Extract the conditional estimates of the fixed effects
701 :     *
702 :     * @param x Pointer to an ssclme object
703 :     *
704 :     * @return a numeric vector containing the conditional estimates of
705 :     * the fixed effects
706 :     */
707 :     SEXP ssclme_fixef(SEXP x)
708 :     {
709 :     SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);
710 :     int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1];
711 :     int j, p = pp1 - 1;
712 :     SEXP val = PROTECT(allocVector(REALSXP, p));
713 :     double
714 :     *RXX = REAL(RXXsl),
715 :     *beta = REAL(val),
716 :     nryyinv; /* negative ryy-inverse */
717 :    
718 :     ssclme_invert(x);
719 :     Memcpy(beta, RXX + p * pp1, p);
720 :     nryyinv = -RXX[pp1*pp1 - 1];
721 :     for (j = 0; j < p; j++) beta[j] /= nryyinv;
722 :     UNPROTECT(1);
723 :     return val;
724 :     }
725 :    
726 :     /**
727 :     * Extract the conditional modes of the random effects.
728 :     *
729 :     * @param x Pointer to an ssclme object
730 :     *
731 :     * @return a vector containing the conditional modes of the random effects
732 :     */
733 :     SEXP ssclme_ranef(SEXP x)
734 :     {
735 : bates 21 SEXP RZXsl = GET_SLOT(x, Matrix_RZXSym),
736 :     GpSl = GET_SLOT(x, Matrix_GpSym);
737 : bates 10 int *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
738 : bates 21 *Gp = INTEGER(GpSl),
739 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
740 :     i, j,
741 : bates 10 n = dims[0],
742 : bates 21 nf = length(GpSl) - 1,
743 : bates 10 pp1 = dims[1],
744 :     p = pp1 - 1;
745 : bates 21 SEXP val = PROTECT(allocVector(VECSXP, nf));
746 : bates 10 double
747 : bates 21 *b = REAL(RZXsl) + n * p,
748 : bates 10 ryyinv; /* ryy-inverse */
749 :    
750 :     ssclme_invert(x);
751 :     ryyinv = REAL(GET_SLOT(x, Matrix_RXXSym))[pp1*pp1 - 1];
752 : bates 21 for (i = 0; i < nf; i++) {
753 : bates 34 int nci = nc[i], Mi = Gp[i+1] - Gp[i];
754 : bates 21 double *mm;
755 :    
756 : bates 34 SET_VECTOR_ELT(val, i, allocMatrix(REALSXP, nci, Mi/nci));
757 : bates 21 mm = Memcpy(REAL(VECTOR_ELT(val, i)), b, Mi);
758 : bates 34 b += Mi;
759 : bates 21 for (j = 0; j < Mi; j++) mm[j] /= ryyinv;
760 :     }
761 : bates 10 UNPROTECT(1);
762 :     return val;
763 :     }
764 : bates 34
765 : bates 10 /**
766 :     * Extract the ML or REML conditional estimate of sigma
767 :     *
768 :     * @param x pointer to an ssclme object
769 :     * @param REML logical scalar - TRUE if REML estimates are requested
770 :     *
771 :     * @return numeric scalar
772 :     */
773 :     SEXP ssclme_sigma(SEXP x, SEXP REML)
774 :     {
775 :     SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);
776 :     int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1],
777 : bates 122 nobs = INTEGER(GET_SLOT(x, Matrix_ncSym))[
778 :     length(GET_SLOT(x, Matrix_GpSym))];
779 : bates 10
780 :     ssclme_invert(x);
781 :     return ScalarReal(1./(REAL(RXXsl)[pp1*pp1 - 1] *
782 :     sqrt((double)(asLogical(REML) ?
783 :     nobs + 1 - pp1 : nobs))));
784 :     }
785 :    
786 :     static
787 :     int coef_length(int nf, const int nc[])
788 :     {
789 :     int i, ans = 0;
790 :     for (i = 0; i < nf; i++) ans += (nc[i] * (nc[i] + 1))/2;
791 :     return ans;
792 :     }
793 :    
794 :     /**
795 :     * Extract the upper triangles of the Omega matrices.
796 : bates 104 * (These aren't "coefficients" but the extractor is
797 : bates 10 * called coef for historical reasons.)
798 :     *
799 :     * @param x pointer to an ssclme object
800 :     *
801 :     * @return numeric vector of the values in the upper triangles of the
802 :     * Omega matrices
803 :     */
804 :     SEXP ssclme_coef(SEXP x)
805 :     {
806 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
807 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
808 :     i, nf = length(Omega), vind;
809 :     SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
810 :     double *vv = REAL(val);
811 :    
812 :     vind = 0;
813 :     for (i = 0; i < nf; i++) {
814 : bates 104 int nci = nc[i];
815 :     if (nci == 1) {
816 :     vv[vind++] = REAL(VECTOR_ELT(Omega, i))[0];
817 :     } else {
818 :     int j, k, odind = vind + nci, ncip1 = nci + 1;
819 :     double *omgi = REAL(VECTOR_ELT(Omega, i));
820 :    
821 :     for (j = 0; j < nci; j++) {
822 :     vv[vind++] = omgi[j * ncip1];
823 :     for (k = j + 1; k < nci; k++) {
824 :     vv[odind++] = omgi[k*nci + j];
825 :     }
826 : bates 10 }
827 : bates 104 vind = odind;
828 : bates 10 }
829 :     }
830 :     UNPROTECT(1);
831 :     return val;
832 :     }
833 :    
834 :     /**
835 : bates 104 * Extract the unconstrained parameters that determine the
836 :     * Omega matrices. (Called coef for historical reasons.)
837 : bates 19 *
838 :     * @param x pointer to an ssclme object
839 :     *
840 : bates 104 * @return numeric vector of unconstrained parameters that determine the
841 : bates 19 * Omega matrices
842 :     */
843 :     SEXP ssclme_coefUnc(SEXP x)
844 :     {
845 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
846 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
847 :     i, nf = length(Omega), vind;
848 :     SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
849 :     double *vv = REAL(val);
850 :    
851 :     vind = 0;
852 :     for (i = 0; i < nf; i++) {
853 :     int nci = nc[i];
854 :     if (nci == 1) {
855 :     vv[vind++] = log(REAL(VECTOR_ELT(Omega, i))[0]);
856 :     } else {
857 :     int j, k, ncip1 = nci + 1, ncisq = nci * nci;
858 :     double *tmp = Memcpy(Calloc(ncisq, double),
859 :     REAL(VECTOR_ELT(Omega, i)), ncisq);
860 :     F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
861 :     if (j) /* should never happen */
862 : bates 104 error("DPOTRF returned error code %d on Omega[[%d]]",
863 :     j, i+1);
864 : bates 19 for (j = 0; j < nci; j++) {
865 :     double diagj = tmp[j * ncip1];
866 :     vv[vind++] = 2. * log(diagj);
867 :     for (k = j + 1; k < nci; k++) {
868 :     tmp[j + k * nci] /= diagj;
869 :     }
870 :     }
871 :     for (j = 0; j < nci; j++) {
872 :     for (k = j + 1; k < nci; k++) {
873 :     vv[vind++] = tmp[j + k * nci];
874 :     }
875 :     }
876 :     Free(tmp);
877 :     }
878 :     }
879 :     UNPROTECT(1);
880 :     return val;
881 :     }
882 :    
883 :     /**
884 : bates 104 * Assign the Omega matrices from the unconstrained parameterization.
885 : bates 19 *
886 :     * @param x pointer to an ssclme object
887 :     * @param coef pointer to an numeric vector of appropriate length
888 :     *
889 :     * @return R_NilValue
890 :     */
891 :     SEXP ssclme_coefGetsUnc(SEXP x, SEXP coef)
892 :     {
893 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
894 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
895 :     cind, i, nf = length(Omega),
896 :     *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
897 :     double *cc = REAL(coef);
898 :    
899 :     if (length(coef) != coef_length(nf, nc) || !isReal(coef))
900 :     error("coef must be a numeric vector of length %d",
901 :     coef_length(nf, nc));
902 :     cind = 0;
903 :     for (i = 0; i < nf; i++) {
904 :     int nci = nc[i];
905 :     if (nci == 1) {
906 :     REAL(VECTOR_ELT(Omega, i))[0] = exp(cc[cind++]);
907 :     } else {
908 :     int odind = cind + nci, /* off-diagonal index */
909 :     j, k,
910 :     ncip1 = nci + 1,
911 :     ncisq = nci * nci;
912 :     double
913 :     *omgi = REAL(VECTOR_ELT(Omega, i)),
914 :     *tmp = Calloc(ncisq, double),
915 : bates 104 diagj, one = 1., zero = 0.;
916 :    
917 : bates 19 memset(omgi, 0, sizeof(double) * ncisq);
918 :     for (j = 0; j < nci; j++) {
919 : bates 104 tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);
920 : bates 19 for (k = j + 1; k < nci; k++) {
921 : bates 104 tmp[k*nci + j] = cc[odind++] * diagj;
922 : bates 19 }
923 :     }
924 : bates 104 F77_CALL(dsyrk)("U", "T", &nci, &nci, &one,
925 :     tmp, &nci, &zero, omgi, &nci);
926 : bates 19 Free(tmp);
927 :     cind = odind;
928 :     }
929 :     }
930 :     status[0] = status[1] = 0;
931 :     return x;
932 :     }
933 :    
934 :     /**
935 : bates 10 * Assign the upper triangles of the Omega matrices.
936 : bates 104 * (Called coef for historical reasons.)
937 : bates 10 *
938 :     * @param x pointer to an ssclme object
939 :     * @param coef pointer to an numeric vector of appropriate length
940 :     *
941 :     * @return R_NilValue
942 :     */
943 :     SEXP ssclme_coefGets(SEXP x, SEXP coef)
944 :     {
945 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
946 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
947 :     cind, i, nf = length(Omega),
948 :     *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
949 :     double *cc = REAL(coef);
950 :    
951 :     if (length(coef) != coef_length(nf, nc) || !isReal(coef))
952 :     error("coef must be a numeric vector of length %d",
953 :     coef_length(nf, nc));
954 :     cind = 0;
955 :     for (i = 0; i < nf; i++) {
956 : bates 104 int nci = nc[i];
957 :     if (nci == 1) {
958 :     REAL(VECTOR_ELT(Omega, i))[0] = cc[cind++];
959 :     } else {
960 :     int j, k, odind = cind + nci, ncip1 = nci + 1;
961 :     double *omgi = REAL(VECTOR_ELT(Omega, i));
962 :    
963 :     for (j = 0; j < nci; j++) {
964 :     omgi[j * ncip1] = cc[cind++];
965 :     for (k = j + 1; k < nci; k++) {
966 :     omgi[k*nci + j] = cc[odind++];
967 :     }
968 : bates 10 }
969 : bates 104 cind = odind;
970 : bates 10 }
971 :     }
972 :     status[0] = status[1] = 0;
973 : bates 11 return x;
974 : bates 10 }
975 :    
976 : bates 11 SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)
977 : bates 10 {
978 :     SEXP
979 :     Omega = GET_SLOT(x, Matrix_OmegaSym),
980 :     RZXsl = GET_SLOT(x, Matrix_RZXSym),
981 :     ncsl = GET_SLOT(x, Matrix_ncSym),
982 :     bVar = GET_SLOT(x, Matrix_bVarSym);
983 :     int
984 :     *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
985 :     *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
986 :     *nc = INTEGER(ncsl),
987 :     *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
988 :     REML = asLogical(REMLp),
989 :     i, info, iter,
990 :     n = dims[0],
991 :     nEM = asInteger(nsteps),
992 :     nf = length(ncsl) - 2,
993 :     nobs = nc[nf + 1],
994 :     p,
995 :     pp1 = dims[1],
996 :     verbose = asLogical(verb);
997 :     double
998 :     *RZX = REAL(RZXsl),
999 :     *dev = REAL(GET_SLOT(x, Matrix_devianceSym)),
1000 :     *b,
1001 :     alpha,
1002 :     one = 1.,
1003 :     zero = 0.;
1004 :    
1005 :     p = pp1 - 1;
1006 :     b = RZX + p * n;
1007 : bates 166 if (verbose) {
1008 :     SEXP coef = PROTECT(ssclme_coef(x));
1009 :     int lc = length(coef); double *cc = REAL(coef);
1010 :    
1011 : bates 168 ssclme_factor(x);
1012 : bates 166 Rprintf(" EM iterations\n");
1013 :     Rprintf("%3d %.3f", 0, dev[REML ? 1 : 0]);
1014 :     for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);
1015 :     Rprintf("\n");
1016 :     UNPROTECT(1);
1017 :     }
1018 :     for (iter = 0; iter < nEM; iter++) {
1019 : bates 10 ssclme_invert(x);
1020 :     for (i = 0; i < nf; i++) {
1021 :     int ki = Gp[i+1] - Gp[i],
1022 :     nci = nc[i],
1023 :     mi = ki/nci;
1024 :     double
1025 :     *vali = REAL(VECTOR_ELT(Omega, i));
1026 :    
1027 :     alpha = ((double)(REML?(nobs-p):nobs))/((double)mi);
1028 :     F77_CALL(dsyrk)("U", "N", &nci, &mi,
1029 :     &alpha, b + Gp[i], &nci,
1030 :     &zero, vali, &nci);
1031 :     alpha = 1./((double) mi);
1032 :     F77_CALL(dsyrk)("U", "N", &nci, &ki,
1033 :     &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,
1034 :     &one, vali, &nci);
1035 :     if (REML) {
1036 : bates 76 int j;
1037 :     for (j = 0; j < p; j++) {
1038 :     F77_CALL(dsyrk)("U", "N", &nci, &mi,
1039 :     &alpha, RZX + Gp[i] + j*n, &nci,
1040 : bates 10 &one, vali, &nci);
1041 : bates 76 }
1042 : bates 10 }
1043 :     F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);
1044 : bates 114 if (info)
1045 :     error("DPOTRF returned error code %d in Omega[[%d]] update",
1046 :     info, i + 1);
1047 : bates 10 F77_CALL(dpotri)("U", &nci, vali, &nci, &info);
1048 : bates 114 if (info)
1049 :     error("DPOTRI returned error code %d in Omega[[%d]] update",
1050 :     info, i + 1);
1051 : bates 10 }
1052 :     status[0] = status[1] = 0;
1053 : bates 166 if (verbose) {
1054 :     SEXP coef = PROTECT(ssclme_coef(x));
1055 :     int lc = length(coef); double *cc = REAL(coef);
1056 : bates 168
1057 :     ssclme_factor(x);
1058 : bates 166 Rprintf("%3d %.3f", iter + 1, dev[REML ? 1 : 0]);
1059 :     for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);
1060 :     Rprintf("\n");
1061 :     UNPROTECT(1);
1062 :     }
1063 : bates 10 }
1064 :     ssclme_factor(x);
1065 :     return R_NilValue;
1066 :     }
1067 : bates 11
1068 : bates 101 SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)
1069 : bates 100 {
1070 :     SEXP
1071 : bates 101 Omega = GET_SLOT(x, Matrix_OmegaSym),
1072 : bates 100 RZXsl = GET_SLOT(x, Matrix_RZXSym),
1073 :     ncsl = GET_SLOT(x, Matrix_ncSym),
1074 :     bVar = GET_SLOT(x, Matrix_bVarSym);
1075 :     int
1076 :     *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1077 :     *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1078 :     *nc = INTEGER(ncsl),
1079 :     REML = asLogical(REMLp),
1080 : bates 104 cind, i, n = dims[0],
1081 : bates 105 nf = length(Omega),
1082 : bates 113 nobs, p, pp1 = dims[1],
1083 : bates 101 uncst = asLogical(Uncp);
1084 : bates 100 double
1085 :     *RZX = REAL(RZXsl),
1086 :     *b,
1087 :     alpha,
1088 :     one = 1.;
1089 : bates 105 SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
1090 : bates 100
1091 : bates 125 ssclme_factor(x);
1092 :     if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0])) {
1093 :     int ncoef = coef_length(nf, nc);
1094 :     for (i = 0; i < ncoef; i++) REAL(ans)[i] = NA_REAL;
1095 :     UNPROTECT(1);
1096 :     return ans;
1097 :     }
1098 : bates 104 nobs = nc[nf + 1];
1099 : bates 100 p = pp1 - 1;
1100 :     b = RZX + p * n;
1101 :     ssclme_invert(x);
1102 : bates 104 cind = 0;
1103 : bates 100 for (i = 0; i < nf; i++) {
1104 : bates 104 int j, ki = Gp[i+1] - Gp[i],
1105 : bates 105 nci = nc[i], ncip1 = nci + 1, ncisq = nci * nci,
1106 : bates 100 mi = ki/nci;
1107 :     double
1108 : bates 105 *chol = Memcpy(Calloc(ncisq, double),
1109 :     REAL(VECTOR_ELT(Omega, i)), ncisq),
1110 :     *tmp = Calloc(ncisq, double);
1111 :    
1112 : bates 100
1113 : bates 105 F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1114 : bates 104 if (j)
1115 :     error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1116 : bates 105 Memcpy(tmp, chol, ncisq);
1117 : bates 104 F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);
1118 :     if (j)
1119 :     error("DPOTRI gave error code %d on Omega[[%d]]", j, i + 1);
1120 : bates 100 alpha = (double) -mi;
1121 :     F77_CALL(dsyrk)("U", "N", &nci, &ki,
1122 :     &one, REAL(VECTOR_ELT(bVar, i)), &nci,
1123 : bates 104 &alpha, tmp, &nci);
1124 : bates 100 alpha = ((double)(REML?(nobs-p):nobs));
1125 :     F77_CALL(dsyrk)("U", "N", &nci, &mi,
1126 :     &alpha, b + Gp[i], &nci,
1127 : bates 104 &one, tmp, &nci);
1128 : bates 100 if (REML) {
1129 :     for (j = 0; j < p; j++) {
1130 :     F77_CALL(dsyrk)("U", "N", &nci, &mi,
1131 :     &one, RZX + Gp[i] + j*n, &nci,
1132 : bates 104 &one, tmp, &nci);
1133 : bates 100 }
1134 :     }
1135 : bates 104 if (nci == 1) {
1136 :     REAL(ans)[cind++] = *tmp *
1137 :     (uncst ? *REAL(VECTOR_ELT(Omega, i)) : 1.);
1138 :     } else {
1139 : bates 106 int k, odind = cind + nci;
1140 : bates 104 if (uncst) {
1141 : bates 106 int ione = 1, kk;
1142 :     double *rr = Calloc(nci, double);
1143 : bates 105 nlme_symmetrize(tmp, nci);
1144 :     for (j = 0; j < nci; j++, cind++) {
1145 : bates 106 for (k = 0; k < nci; k++) rr[k] = chol[j + k*nci];
1146 :     REAL(ans)[cind] = 0.;
1147 : bates 105 for (k = j; k < nci; k++) {
1148 :     for (kk = j; kk < nci; kk++) {
1149 : bates 106 REAL(ans)[cind] += rr[k] * rr[kk] *
1150 : bates 105 tmp[kk * nci + k];
1151 :     }
1152 :     }
1153 : bates 106 for (k = 0; k < nci; k++) rr[k] *= rr[j];
1154 : bates 105 for (k = j + 1; k < nci; k++) {
1155 : bates 106 REAL(ans)[odind++] =
1156 :     F77_CALL(ddot)(&nci, rr, &ione, tmp + k, &nci) +
1157 :     F77_CALL(ddot)(&nci, rr, &ione,
1158 :     tmp + k*nci, &ione);
1159 : bates 105 }
1160 :     }
1161 : bates 106 Free(rr);
1162 : bates 101 } else {
1163 : bates 104 for (j = 0; j < nci; j++) {
1164 :     REAL(ans)[cind++] = tmp[j * ncip1];
1165 :     for (k = j + 1; k < nci; k++) {
1166 :     REAL(ans)[odind++] = tmp[k*nci + j] * 2.;
1167 :     }
1168 :     }
1169 : bates 101 }
1170 : bates 105 cind = odind;
1171 : bates 101 }
1172 : bates 105 Free(tmp); Free(chol);
1173 : bates 100 }
1174 :     UNPROTECT(1);
1175 :     return ans;
1176 :     }
1177 :    
1178 : bates 164 SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats, SEXP useRf)
1179 : bates 28 {
1180 :     SEXP val, b;
1181 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1182 :     i, ione = 1, nf = length(facs), nobs, p;
1183 :     double *vv, one = 1.0, zero = 0.0;
1184 :    
1185 :     if (nf < 1)
1186 :     error("null factor list passed to ssclme_fitted");
1187 :     nobs = length(VECTOR_ELT(facs, 0));
1188 :     val = PROTECT(allocVector(REALSXP, nobs));
1189 :     vv = REAL(val);
1190 :     p = nc[nf] - 1;
1191 :     if (p > 0) {
1192 :     F77_CALL(dgemm)("N", "N", &nobs, &ione, &p, &one,
1193 :     REAL(VECTOR_ELT(mmats, nf)), &nobs,
1194 :     REAL(PROTECT(ssclme_fixef(x))), &p,
1195 :     &zero, vv, &nobs);
1196 :     UNPROTECT(1);
1197 :     } else {
1198 :     memset(vv, 0, sizeof(double) * nobs);
1199 :     }
1200 : bates 164 if (asLogical(useRf)) {
1201 :     b = PROTECT(ssclme_ranef(x));
1202 :     for (i = 0; i < nf; i++) {
1203 :     int *ff = INTEGER(VECTOR_ELT(facs, i)), j, nci = nc[i];
1204 :     double *bb = REAL(VECTOR_ELT(b, i)),
1205 :     *mm = REAL(VECTOR_ELT(mmats, i));
1206 :     for (j = 0; j < nobs; ) {
1207 :     int nn = 1, lev = ff[j];
1208 :     /* check for adjacent rows with same factor level */
1209 :     while (ff[j + nn] == lev) nn++;
1210 :     F77_CALL(dgemm)("N", "N", &nn, &ione, &nci,
1211 :     &one, mm + j, &nobs,
1212 :     bb + (lev - 1) * nci, &nci,
1213 :     &one, vv + j, &nobs);
1214 :     j += nn;
1215 :     }
1216 : bates 28 }
1217 : bates 164 UNPROTECT(1);
1218 : bates 28 }
1219 : bates 164 UNPROTECT(1);
1220 : bates 28 return val;
1221 :     }
1222 : bates 45
1223 : bates 122 /**
1224 :     * Return the unscaled variances
1225 :     *
1226 :     * @param x pointer to an ssclme object
1227 :     *
1228 :     * @return a list similar to the Omega list with the unscaled variances
1229 :     */
1230 :     SEXP ssclme_variances(SEXP x)
1231 : bates 45 {
1232 : bates 122 SEXP Omg = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym)));
1233 : bates 45 int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1234 : bates 57 i, nf = length(Omg);
1235 : bates 45
1236 :     for (i = 0; i < nf; i++) {
1237 : bates 57 double *mm = REAL(VECTOR_ELT(Omg, i));
1238 : bates 122 int j, nci = nc[i];
1239 : bates 45
1240 :     F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);
1241 :     if (j) /* shouldn't happen */
1242 :     error("DPOTRF returned error code %d on Omega[%d]",
1243 :     j, i + 1);
1244 :     F77_CALL(dpotri)("U", &nci, mm, &nci, &j);
1245 :     if (j) /* shouldn't happen */
1246 :     error("DTRTRI returned error code %d on Omega[%d]",
1247 :     j, i + 1);
1248 : bates 122 nlme_symmetrize(mm, nci);
1249 : bates 45 }
1250 : bates 122 UNPROTECT(1);
1251 :     return Omg;
1252 : bates 45 }
1253 : bates 100
1254 : bates 164 #define slot_dup(sym) SET_SLOT(ans, sym, duplicate(GET_SLOT(x, sym)))
1255 :    
1256 : bates 114 SEXP ssclme_collapse(SEXP x)
1257 :     {
1258 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("ssclme"))),
1259 :     Omega = GET_SLOT(x, Matrix_OmegaSym),
1260 :     Dim = GET_SLOT(x, Matrix_DimSym);
1261 : bates 164 int nf = length(Omega), nz = INTEGER(Dim)[1];
1262 : bates 114
1263 : bates 164 slot_dup(Matrix_DSym);
1264 :     slot_dup(Matrix_DIsqrtSym);
1265 :     slot_dup(Matrix_DimSym);
1266 :     slot_dup(Matrix_GpSym);
1267 :     slot_dup(Matrix_LiSym);
1268 :     slot_dup(Matrix_LpSym);
1269 :     slot_dup(Matrix_LxSym);
1270 :     slot_dup(Matrix_OmegaSym);
1271 :     slot_dup(Matrix_ParentSym);
1272 :     slot_dup(Matrix_bVarSym);
1273 :     slot_dup(Matrix_devianceSym);
1274 :     slot_dup(Matrix_devCompSym);
1275 :     slot_dup(Matrix_iSym);
1276 :     slot_dup(Matrix_ncSym);
1277 :     slot_dup(Matrix_statusSym);
1278 :     slot_dup(Matrix_pSym);
1279 :     slot_dup(Matrix_xSym);
1280 : bates 114 INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;
1281 :     SET_SLOT(ans, Matrix_XtXSym, allocMatrix(REALSXP, 1, 1));
1282 : bates 122 REAL(GET_SLOT(ans, Matrix_XtXSym))[0] = NA_REAL;
1283 : bates 114 SET_SLOT(ans, Matrix_RXXSym, allocMatrix(REALSXP, 1, 1));
1284 : bates 122 REAL(GET_SLOT(ans, Matrix_RXXSym))[0] = NA_REAL;
1285 : bates 114 SET_SLOT(ans, Matrix_ZtXSym, allocMatrix(REALSXP, nz, 1));
1286 :     SET_SLOT(ans, Matrix_RZXSym, allocMatrix(REALSXP, nz, 1));
1287 : bates 115 LOGICAL(GET_SLOT(ans, Matrix_statusSym))[0] = 0;
1288 : bates 114 UNPROTECT(1);
1289 :     return ans;
1290 :     }
1291 :    
1292 : bates 164 SEXP ssclme_to_lme(SEXP call, SEXP facs, SEXP x, SEXP model, SEXP REML,
1293 :     SEXP rep, SEXP fitted, SEXP residuals)
1294 :     {
1295 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lme")));
1296 :    
1297 :     SET_SLOT(ans, install("call"), call);
1298 :     SET_SLOT(ans, install("facs"), facs);
1299 :     SET_SLOT(ans, Matrix_xSym, x);
1300 :     SET_SLOT(ans, install("model"), model);
1301 :     SET_SLOT(ans, install("REML"), REML);
1302 :     SET_SLOT(ans, install("rep"), rep);
1303 :     SET_SLOT(ans, install("fitted"), fitted);
1304 :     SET_SLOT(ans, install("residuals"), residuals);
1305 :     UNPROTECT(1);
1306 :     return ans;
1307 :     }

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