SCM Repository

[matrix] Diff of /pkg/src/ssclme.c
 [matrix] / pkg / src / ssclme.c

Diff of /pkg/src/ssclme.c

revision 22, Sat Mar 27 18:17:00 2004 UTC revision 176, Fri May 28 13:23:44 2004 UTC
# Line 1  Line 1
1  #include "ssclme.h"  #include "ssclme.h"
2
3  /**  #define slot_dup(dest, src, sym)  SET_SLOT(dest, sym, duplicate(GET_SLOT(src, sym)))
* Check for a nested series of grouping factors in the sparse,
*  symmetric representation of the pairwise cross-tabulations.
*
* @param n size of pairwise cross-tabulation matrix
* @param nf number of groups of columns in pairwise cross-tabulation
* @param upper non-zero if the upper triangle is stored
* @param Ap array of pointers to columns
* @param Ai row indices
* @param Gp array of pointers to groups
*
* @return 0 for non-nested groups, 1 for nested groups
*/
static
int ctab_isNested(int n, int nf, int upper,
const int Ap[], const int Ai[], const int Gp[])
{
if (nf > 1) {  /* single factor always nested */
int  i;
if (upper) {
int *nnz = (int *) R_alloc(n, sizeof(int)), nz = Ap[n];
/* count number of nonzeros in each row */
for (i = 0; i < n; i++) nnz[i] = 0;
for (i = 0; i < nz; i++) nnz[Ai[i]]++;
for (i = 0; i < nf; i++) {
int j, p2 = Gp[i+1], target = nf - i;
for (j = Gp[i]; j < p2; j++) {
if (nnz[j] != target) return 0;
}
}
} else {                /* lower triangle - the easy case */
for (i = 0; i < nf; i++) {
int j, p2 = Gp[i+1], target = nf - i;
for (j = Gp[i]; j < p2; j++) {
if ((Ap[j+1] - Ap[j]) != target)
return 0;
}
}
}
}
return 1;
}
4
5  /**  /**
6   * Determine if a fill-reducing permutation is needed for the pairwise   * Using the sscCrosstab object from the grouping factors, generate
7   * cross-tabulation matrix.  If so, determine such a permutation   * the slots in an ssclme object related to the symmetric sparse
8   * (using Metis) then separate the groups.   * matrix representation of Z'Z.  If the model matrices for the
9     * grouping factors have only one column each then the structure can
10     * be copied, otherwise it must be generated from the sscCrosstab and
11     * the number of columns per grouping factor.
12   *   *
13   * @param ctab pointer to a pairwise cross-tabulation object   * @param nf number of factors
14   *   * @param nc vector of length nf+2 with number of columns in model matrices
15   * @return pointer to an integer R vector.   * @param ctab pointer to the sscCrosstab object
16     * @param ssc pointer to an ssclme object to be filled out
17   */   */
18  static  static
SEXP ctab_permute(SEXP ctab)
{
SEXP val, GpSl = GET_SLOT(ctab, Matrix_GpSym);
int *Ai = INTEGER(GET_SLOT(ctab, Matrix_iSym)),
*Ap = INTEGER(GET_SLOT(ctab, Matrix_pSym)),
*Gp = INTEGER(GpSl),
*perm,
*work,
i,
j,
n = INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],
nf = length(GpSl) - 1,
nz = Ap[n],             /* number of non-zeros */
pos;

if (ctab_isNested(n, nf, 1, Ap, Ai, Gp))
return allocVector(INTSXP, 0);
val =  allocVector(INTSXP, n);
perm = INTEGER(val);
work = (int *) R_alloc(n, sizeof(int));
ssc_metis_order(n, nz, Ap, Ai, work, perm); /* perm gets inverse perm */
/* work now contains desired permutation but with groups scrambled */

/* copy work into perm preserving the order of the groups */
pos = 0;            /* position in new permutation */
for (i = 0; i < nf; i++) {
for (j = 0; j < n; j++) {
int jj = work[j];
if (Gp[i] <= jj && jj < Gp[i+1]) {
perm[pos] = jj;
pos++;
}
}
}
return val;
}

static
19  void ssclme_copy_ctab(int nf, const int nc[], SEXP ctab, SEXP ssc)  void ssclme_copy_ctab(int nf, const int nc[], SEXP ctab, SEXP ssc)
20  {  {
21      int *snc, i, copyonly = 1;      int *snc, i, copyonly = 1;
22
23      for (i = 0; i < nf; i++) {      SET_SLOT(ssc, Matrix_ncSym, allocVector(INTSXP, nf + 2));
24          if (nc[i] > 1) copyonly = 0;      snc = INTEGER(GET_SLOT(ssc, Matrix_ncSym));
25        for (i = 0; i <= nf; i++) {
26            snc[i] = nc[i];
27            if (nc[i] > 1 && i < nf) copyonly = 0;
28      }      }
29      if (copyonly) {      if (copyonly) {
30          SET_SLOT(ssc, Matrix_pSym, duplicate(GET_SLOT(ctab, Matrix_pSym)));          slot_dup(ssc, ctab, Matrix_pSym);
31          SET_SLOT(ssc, Matrix_iSym, duplicate(GET_SLOT(ctab, Matrix_iSym)));          slot_dup(ssc, ctab, Matrix_iSym);
32          SET_SLOT(ssc, Matrix_xSym, duplicate(GET_SLOT(ctab, Matrix_xSym)));          slot_dup(ssc, ctab, Matrix_xSym);
33          SET_SLOT(ssc, Matrix_DimSym,          slot_dup(ssc, ctab, Matrix_DimSym);
34                   duplicate(GET_SLOT(ctab, Matrix_DimSym)));          slot_dup(ssc, ctab, Matrix_GpSym);
35          SET_SLOT(ssc, Matrix_GpSym, duplicate(GET_SLOT(ctab, Matrix_GpSym)));          return;
36      } else {      }
37        {
38          int          int
39              *GpIn = INTEGER(GET_SLOT(ctab, Matrix_GpSym)),              *GpIn = INTEGER(GET_SLOT(ctab, Matrix_GpSym)),
40              *GpOut,              *GpOut,
# Line 131  Line 60
60              }              }
61          }          }
62          nOut = GpOut[nf];       /* size of output matrix */          nOut = GpOut[nf];       /* size of output matrix */
63          SET_SLOT(ssc, Matrix_DimSym,          SET_SLOT(ssc, Matrix_DimSym, allocVector(INTSXP, 2));
duplicate(GET_SLOT(ctab, Matrix_DimSym)));
64          dims = INTEGER(GET_SLOT(ssc, Matrix_DimSym));          dims = INTEGER(GET_SLOT(ssc, Matrix_DimSym));
65          dims[0] = dims[1] = nOut;          dims[0] = dims[1] = nOut;
66          SET_SLOT(ssc, Matrix_pSym, allocVector(INTSXP, nOut + 1));          SET_SLOT(ssc, Matrix_pSym, allocVector(INTSXP, nOut + 1));
# Line 163  Line 91
91          for (i = 0; i < nf; i++) { /* fill in the rows */          for (i = 0; i < nf; i++) { /* fill in the rows */
92              int j, jj, nci = nc[i], p2 = GpIn[i+1];              int j, jj, nci = nc[i], p2 = GpIn[i+1];
93              for (j = GpIn[i]; j < p2; j++) { /* first col for input col */              for (j = GpIn[i]; j < p2; j++) { /* first col for input col */
94                  int ii = AiIn[j], mj = map[j], ncci = ncc[ii],                  int k, mj = map[j], p3 = ApIn[j+1], pos = ApOut[mj];
95                      pos = ApOut[mj];                  for (k = ApIn[j]; k < p3; k++) {
96                        int ii = AiIn[k], ncci = ncc[ii];
97                  AiOut[pos++] = map[ii];                  AiOut[pos++] = map[ii];
98                  if (ii < j) {   /* above the diagonal */                  if (ii < j) {   /* above the diagonal */
99                      for (jj = 1; jj < ncci; jj++) {                      for (jj = 1; jj < ncci; jj++) {
# Line 176  Line 105
105                                                  * another diagonal element */                                                  * another diagonal element */
106                      int mjj = mj + jj, pj = ApOut[mjj], pjm1 = ApOut[mjj-1];                      int mjj = mj + jj, pj = ApOut[mjj], pjm1 = ApOut[mjj-1];
107                      Memcpy(AiOut + pj, AiOut + pjm1, pj - pjm1);                      Memcpy(AiOut + pj, AiOut + pjm1, pj - pjm1);
108                      AiOut[ApOut[mjj + 1] - 1] = mjj; /* maybe mjj-1? */                          AiOut[ApOut[mjj + 1] - 1] = mjj;
109                  }                  }
110              }              }
111          }          }
Free(map); Free(ncc);
112      }      }
113      SET_SLOT(ssc, Matrix_ncSym, allocVector(INTSXP, nf + 2));          Free(map); Free(ncc);
snc = INTEGER(GET_SLOT(ssc, Matrix_ncSym));
for (i = 0; i <= nf; i++) {
snc[i] = nc[i];
114      }      }
115  }  }
116
117  static  /**
118  void ssclme_fill_LIp(int n, const int Parent[], int LIp[])   * Calculate and store the maximum number of off-diagonal elements in
119     * the inverse of L, based on the elimination tree.  The maximum is
120     * itself stored in the Parent array.  (FIXME: come up with a better design.)
121     *
122     * @param n number of columns in the matrix
123     * @param Parent elimination tree for the matrix
124     */
125    static void ssclme_calc_maxod(int n, int Parent[])
126  {  {
127      int *sz = Calloc(n, int), i;      int *sz = Calloc(n, int), i, mm = -1;
128      for (i = n - 1; i >= 0; i--) {      for (i = n - 1; i >= 0; i--) {
129          sz[i] = (Parent[i] < 0) ? 0 : 1 + sz[Parent[i]];          sz[i] = (Parent[i] < 0) ? 0 : (1 + sz[Parent[i]]);
130            if (sz[i] > mm) mm = sz[i];
131      }      }
132      LIp[0] = 0;      Parent[n] = mm;
for (i = 0; i < n; i++) LIp[i+1] = LIp[i] + sz[i];
133      Free(sz);      Free(sz);
134  }  }
135
136  static  /**
137  void ssclme_fill_LIi(int n, const int Parent[], const int LIp[], int LIi[])   * Create an ssclme object from a list of grouping factors, sorted in
138  {   * order of non-increasing numbers of levels, and an integer vector of
139      int i;   * the number of columns in the model matrices.  There is one more
140      for (i = n; i > 0; i--) {   * element in ncv than in facs.  The last element is the number of
141          int im1 = i - 1, Par = Parent[im1];   * columns in the model matrix for the fixed effects plus the
142          if (Par >= 0) {   * response.  (i.e. p+1)
143              LIi[LIp[im1]] = Par;   *
144              Memcpy(LIi + LIp[im1] + 1, LIi + LIp[Par],   * @param facs pointer to a list of grouping factors
145                     LIp[Par + 1] - LIp[Par]);   * @param ncv pointer to an integer vector of number of columns per model matrix
146          }   *
147      }   * @return pointer to an ssclme object
148  }   */

149  SEXP  SEXP
150  ssclme_create(SEXP facs, SEXP ncv, SEXP threshold)  ssclme_create(SEXP facs, SEXP ncv)
151  {  {
152      SEXP ctab, nms, ssc, tmp,      SEXP ctab, nms, ssc, tmp,
153          val = PROTECT(allocVector(VECSXP, 2)),          val = PROTECT(allocVector(VECSXP, 2)),
154          dd = PROTECT(allocVector(INTSXP, 3));   /* dimensions of 3-D arrays */          dd = PROTECT(allocVector(INTSXP, 3));   /* dimensions of 3-D arrays */
155      int *Ai, *Ap, *Gp, *LIp, *Lp, *Parent,      int *Ai, *Ap, *Gp, *Lp, *Parent,
156          *nc, Lnz, i, nf = length(facs), nzcol, pp1,          *nc, Lnz, i, nf = length(facs), nzcol, pp1,
157          *dims = INTEGER(dd);          *dims = INTEGER(dd);
158
# Line 232  Line 163
163      ssc = VECTOR_ELT(val, 0);      ssc = VECTOR_ELT(val, 0);
164                                  /* Pairwise cross-tabulation */                                  /* Pairwise cross-tabulation */
165      ctab = PROTECT(sscCrosstab(facs, ScalarLogical(1)));      ctab = PROTECT(sscCrosstab(facs, ScalarLogical(1)));
166      SET_VECTOR_ELT(val, 1, ctab_permute(ctab));      SET_VECTOR_ELT(val, 1, sscCrosstab_groupedPerm(ctab));
167      if (length(VECTOR_ELT(val, 1)) > 0) {/* Fill-reducing permutation */      if (length(VECTOR_ELT(val, 1)) > 0) {/* Fill-reducing permutation */
168          ssc_symbolic_permute(INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],          ssc_symbolic_permute(INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],
169                               1, INTEGER(VECTOR_ELT(val, 1)),                               1, INTEGER(VECTOR_ELT(val, 1)),
# Line 254  Line 185
185      SET_SLOT(ssc, Matrix_RXXSym, allocMatrix(REALSXP, pp1, pp1));      SET_SLOT(ssc, Matrix_RXXSym, allocMatrix(REALSXP, pp1, pp1));
186      SET_SLOT(ssc, Matrix_ZtXSym, allocMatrix(REALSXP, nzcol, pp1));      SET_SLOT(ssc, Matrix_ZtXSym, allocMatrix(REALSXP, nzcol, pp1));
187      SET_SLOT(ssc, Matrix_RZXSym, allocMatrix(REALSXP, nzcol, pp1));      SET_SLOT(ssc, Matrix_RZXSym, allocMatrix(REALSXP, nzcol, pp1));
188          /* Zero the symmetric matrices (for cosmetic reasons only). */                                  /* Zero symmetric matrices (cosmetic) */
189      memset(REAL(GET_SLOT(ssc, Matrix_XtXSym)), 0,      memset(REAL(GET_SLOT(ssc, Matrix_XtXSym)), 0,
190             sizeof(double) * pp1 * pp1);             sizeof(double) * pp1 * pp1);
191      memset(REAL(GET_SLOT(ssc, Matrix_RXXSym)), 0,      memset(REAL(GET_SLOT(ssc, Matrix_RXXSym)), 0,
192             sizeof(double) * pp1 * pp1);             sizeof(double) * pp1 * pp1);
193      SET_SLOT(ssc, Matrix_LpSym, allocVector(INTSXP, nzcol + 1));      SET_SLOT(ssc, Matrix_LpSym, allocVector(INTSXP, nzcol + 1));
194      Lp = INTEGER(GET_SLOT(ssc, Matrix_LpSym));      Lp = INTEGER(GET_SLOT(ssc, Matrix_LpSym));
195      SET_SLOT(ssc, Matrix_ParentSym, allocVector(INTSXP, nzcol));      SET_SLOT(ssc, Matrix_ParentSym, allocVector(INTSXP, nzcol + 1));
196      Parent = INTEGER(GET_SLOT(ssc, Matrix_ParentSym));      Parent = INTEGER(GET_SLOT(ssc, Matrix_ParentSym));
197      SET_SLOT(ssc, Matrix_DSym, allocVector(REALSXP, nzcol));      SET_SLOT(ssc, Matrix_DSym, allocVector(REALSXP, nzcol));
198      SET_SLOT(ssc, Matrix_DIsqrtSym, allocVector(REALSXP, nzcol));      SET_SLOT(ssc, Matrix_DIsqrtSym, allocVector(REALSXP, nzcol));
# Line 269  Line 200
200                   (int *) R_alloc(nzcol, sizeof(int)), /* Lnz */                   (int *) R_alloc(nzcol, sizeof(int)), /* Lnz */
201                   (int *) R_alloc(nzcol, sizeof(int)), /* Flag */                   (int *) R_alloc(nzcol, sizeof(int)), /* Flag */
202                   (int *) NULL, (int *) NULL); /* P & Pinv */                   (int *) NULL, (int *) NULL); /* P & Pinv */
203        ssclme_calc_maxod(nzcol, Parent);
204      Lnz = Lp[nzcol];      Lnz = Lp[nzcol];
205      SET_SLOT(ssc, Matrix_LiSym, allocVector(INTSXP, Lnz));      SET_SLOT(ssc, Matrix_LiSym, allocVector(INTSXP, Lnz));
206      SET_SLOT(ssc, Matrix_LxSym, allocVector(REALSXP, Lnz));      SET_SLOT(ssc, Matrix_LxSym, allocVector(REALSXP, Lnz));
# Line 306  Line 238
238          memset(REAL(VECTOR_ELT(tmp, i)), 0,          memset(REAL(VECTOR_ELT(tmp, i)), 0,
239                 sizeof(double) * nci * nci * mi);                 sizeof(double) * nci * nci * mi);
240      }      }
SET_SLOT(ssc, Matrix_LIpSym, allocVector(INTSXP, nzcol + 1));
LIp = INTEGER(GET_SLOT(ssc, Matrix_LIpSym));
ssclme_fill_LIp(nzcol, Parent, LIp);
if (asInteger(threshold) > (Lnz = LIp[nzcol])) {
SET_SLOT(ssc, Matrix_LIiSym, allocVector(INTSXP, Lnz));
ssclme_fill_LIi(nzcol, Parent, LIp,
INTEGER(GET_SLOT(ssc, Matrix_LIiSym)));
SET_SLOT(ssc, Matrix_LIxSym, allocVector(REALSXP, Lnz));
memset(REAL(GET_SLOT(ssc, Matrix_LIxSym)), 0,
sizeof(double) * Lnz);
}
241      UNPROTECT(2);      UNPROTECT(2);
242      return val;      return val;
243  }  }
244
245    /**
246     * Copy information on Z'Z accumulated in the bVar array to Z'Z
247     *
248     * @param ncj number of columns in this block
249     * @param Gpj initial column for this group
250     * @param Gpjp initial column for the next group
251     * @param bVj pointer to the ncj x ncj x mj array to be filled
252     * @param Ap column pointer array for Z'Z
253     * @param Ai row indices for Z'Z
254     * @param Ax elements of Z'Z
255     */
256  static  static
257  void bVj_to_A(int ncj, int Gpj, int Gpjp, const double bVj[],  void bVj_to_A(int ncj, int Gpj, int Gpjp, const double bVj[],
258                const int Ap[], const int Ai[], double Ax[])                const int Ap[], const int Ai[], double Ax[])
# Line 337  Line 269
269      }      }
270  }  }
271
272    /**
273     * Copy the dimnames from the list of grouping factors and the model
274     * matrices for the grouping factors into the appropriate parts of the
275     * ssclme object.
276     *
277     * @param x pointer to an ssclme object
278     * @param facs pointer to a list of factors
279     * @param mmats pointer to a list of model matrices
280     *
281     * @return NULL
282     */
283  SEXP  SEXP
284  ssclme_transfer_dimnames(SEXP x, SEXP facs, SEXP mmats)  ssclme_transfer_dimnames(SEXP x, SEXP facs, SEXP mmats)
285  {  {
# Line 362  Line 305
305      return R_NilValue;      return R_NilValue;
306  }  }
307
308    /**
309     * Update the numerical entries x, ZtX, and XtX in an ssclme object
310     * according to a set of model matrices.
311     *
312     * @param x pointer to an ssclme object
313     * @param facs pointer to a list of grouping factors
314     * @param mmats pointer to a list of model matrices
315     *
316     * @return NULL
317     */
318  SEXP  SEXP
319  ssclme_update_mm(SEXP x, SEXP facs, SEXP mmats)  ssclme_update_mm(SEXP x, SEXP facs, SEXP mmats)
320  {  {
# Line 371  Line 324
324          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
325          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
326          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
327            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
328          i, j, k,          i, j, k,
329          ione = 1,          ione = 1,
330          nf = length(mmats) - 1,          nf = length(mmats) - 1,
# Line 434  Line 388
388          for (k = j+1; k < nf; k++) { /* off-diagonals */          for (k = j+1; k < nf; k++) { /* off-diagonals */
389              int *fpk = INTEGER(VECTOR_ELT(facs, k)),              int *fpk = INTEGER(VECTOR_ELT(facs, k)),
390                  *Apk = Ap + Gp[k],                  *Apk = Ap + Gp[k],
391                  nck = nc[k];                  nck = nc[k],
392                    scalar = ncj == 1 && nck == 1;
393              double              double
394                  *Zk = Z[k];                  *Zk = Z[k], *work;
395                if (!scalar) work = Calloc(ncj * nck, double);
396              for (i = 0; i < nobs; i++) {              for (i = 0; i < nobs; i++) {
397                  int ii, ind = -1, fpji = fpj[i] - 1,                  int ii, ind = -1, fpji = fpj[i] - 1,
398                      row = Gp[j] + fpji * ncj,                      row = Gp[j] + fpji * ncj,
399                      fpki = fpk[i] - 1,                      fpki = fpk[i] - 1,
400                      lastind = Apk[fpki + 1];                      lastind = Apk[fpki*nck + 1];
401                  for (ii = Apk[fpki]; ii < lastind; ii++) {                  for (ii = Apk[fpki*nck]; ii < lastind; ii++) {
402                      if (Ai[ii] == row) {                      if (Ai[ii] == row) {
403                          ind = ii;                          ind = ii;
404                          break;                          break;
405                      }                      }
406                  }                  }
407                  if (ind < 0) error("logic error in ssclme_update_mm");                  if (ind < 0) error("logic error in ssclme_update_mm");
408                  if (Ncj || nck > 1) {                  if (scalar) {   /* update scalars directly */
409                                  /* FIXME: run a loop to update */                      Ax[ind] += Zj[i] * Zk[i];
410                      error("code not yet written");                  } else {
411                  } else {        /* update scalars directly */                      int jj, offset = ind - Apk[fpki * nck];
412                      Ax[ind] += Zj[fpji] * Zk[fpki];                      F77_CALL(dgemm)("T", "N", &ncj, &nck, &ione, &one,
413                                        Zj + i, &nobs, Zk + i, &nobs,
414                                        &zero, work, &ncj);
415                        for (jj = 0; jj < nck; jj++) {
416                            ind = Apk[fpki * nck + jj] + offset;
417                            if (Ai[ind] != row)
418                                error("logic error in ssclme_update_mm");
419                            for (ii = 0; ii < ncj; ii++) {
420                                Ax[ind++] += work[jj * ncj + ii];
421                  }                  }
422              }              }
423          }          }
424      }      }
425                if (!scalar) Free(work);
426            }
427        }
428      Free(Z);      Free(Z);
429      ssclme_transfer_dimnames(x, facs, mmats);      ssclme_transfer_dimnames(x, facs, mmats);
430        status[0] = status[1] = 0;
431      return R_NilValue;      return R_NilValue;
432  }  }
433
434  SEXP ssclme_inflate_and_factor(SEXP lme)  /**
435     * Inflate Z'Z according to Omega and create the factorization LDL'
436     *
437     * @param x pointer to an ssclme object
438     *
439     * @return NULL
440     */
441    SEXP ssclme_inflate_and_factor(SEXP x)
442  {  {
443      SEXP      SEXP
444          GpSlot = GET_SLOT(lme, Matrix_GpSym),          GpSlot = GET_SLOT(x, Matrix_GpSym),
445          Omega = GET_SLOT(lme, Matrix_OmegaSym);          Omega = GET_SLOT(x, Matrix_OmegaSym);
446      int n = INTEGER(GET_SLOT(lme, Matrix_DimSym))[1];      int n = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];
447      int      int
448          *Ai = INTEGER(GET_SLOT(lme, Matrix_iSym)),          *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
449          *Ap = INTEGER(GET_SLOT(lme, Matrix_pSym)),          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
450          *Flag = Calloc(n, int),          *Flag = Calloc(n, int),
451          *Gp = INTEGER(GpSlot),          *Gp = INTEGER(GpSlot),
452          *Lnz = Calloc(n, int),          *Lnz = Calloc(n, int),
453          *Pattern = Calloc(n, int),          *Pattern = Calloc(n, int),
454          *nc = INTEGER(GET_SLOT(lme, Matrix_ncSym)),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
455          j,          j,
456          nf = length(GpSlot) - 1;          nf = length(GpSlot) - 1;
457      double      double
458          *D = REAL(GET_SLOT(lme, Matrix_DSym)),          *D = REAL(GET_SLOT(x, Matrix_DSym)),
459          *DIsqrt = REAL(GET_SLOT(lme, Matrix_DIsqrtSym)),          *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
460          *Y = Calloc(n, double),          *Y = Calloc(n, double),
461          *xcp = Calloc(Ap[n], double);          *xcp = Calloc(Ap[n], double);
462
463      Memcpy(xcp, REAL(GET_SLOT(lme, Matrix_xSym)), Ap[n]);      Memcpy(xcp, REAL(GET_SLOT(x, Matrix_xSym)), Ap[n]);
464      for (j = 0; j < nf; j++) {      for (j = 0; j < nf; j++) {
465          int  diag, i, ii, k, G2 = Gp[j + 1], ncj = nc[j];          int  diag, i, ii, k, G2 = Gp[j + 1], ncj = nc[j];
466          double *omgj = REAL(VECTOR_ELT(Omega, j));          double *omgj = REAL(VECTOR_ELT(Omega, j));
# Line 504  Line 478
478          }          }
479      }      }
480      j = ldl_numeric(n, Ap, Ai, xcp,      j = ldl_numeric(n, Ap, Ai, xcp,
481                      INTEGER(GET_SLOT(lme, Matrix_LpSym)),                      INTEGER(GET_SLOT(x, Matrix_LpSym)),
482                      INTEGER(GET_SLOT(lme, Matrix_ParentSym)),                      INTEGER(GET_SLOT(x, Matrix_ParentSym)),
483                      Lnz, INTEGER(GET_SLOT(lme, Matrix_LiSym)),                      Lnz, INTEGER(GET_SLOT(x, Matrix_LiSym)),
484                      REAL(GET_SLOT(lme, Matrix_LxSym)),                      REAL(GET_SLOT(x, Matrix_LxSym)),
485                      D, Y, Pattern, Flag,                      D, Y, Pattern, Flag,
486                      (int *) NULL, (int *) NULL); /* P & Pinv */                      (int *) NULL, (int *) NULL); /* P & Pinv */
487      if (j != n)      if (j != n)
# Line 518  Line 492
492      return R_NilValue;      return R_NilValue;
493  }  }
494
495  SEXP ssclme_factor(SEXP lme)
496    /**
497     * If status[["factored"]] is FALSE, create and factor Z'Z+Omega, then
498     * create RZX and RXX, the deviance components, and the value of the
499     * deviance for both ML and REML.
500     *
501     * @param x pointer to an ssclme object
502     *
503     * @return NULL
504     */
505    SEXP ssclme_factor(SEXP x)
506  {  {
507      int *status = LOGICAL(GET_SLOT(lme, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
508
509      if (!status[0]) {      if (!status[0]) {
510          SEXP          SEXP
511              GpSlot = GET_SLOT(lme, Matrix_GpSym),              GpSlot = GET_SLOT(x, Matrix_GpSym),
512              Omega = GET_SLOT(lme, Matrix_OmegaSym);              Omega = GET_SLOT(x, Matrix_OmegaSym);
513          int          int
514              *Gp = INTEGER(GpSlot),              *Gp = INTEGER(GpSlot),
515              *Li = INTEGER(GET_SLOT(lme, Matrix_LiSym)),              *Li = INTEGER(GET_SLOT(x, Matrix_LiSym)),
516              *Lp = INTEGER(GET_SLOT(lme, Matrix_LpSym)),              *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
517              *nc = INTEGER(GET_SLOT(lme, Matrix_ncSym)),              *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
518              i,              i,
519              n = INTEGER(GET_SLOT(lme, Matrix_DimSym))[1],              n = INTEGER(GET_SLOT(x, Matrix_DimSym))[1],
520              nf = length(GpSlot) - 1,              nf = length(GpSlot) - 1,
521              nobs = nc[nf + 1],              nobs = nc[nf + 1],
522              nreml = nobs + 1 - nc[nf],              nreml = nobs + 1 - nc[nf],
523              pp1 = nc[nf],              pp1 = nc[nf],
524              pp2 = pp1 + 1;              pp2 = pp1 + 1;
525          double          double
526              *D = REAL(GET_SLOT(lme, Matrix_DSym)),              *D = REAL(GET_SLOT(x, Matrix_DSym)),
527              *DIsqrt = REAL(GET_SLOT(lme, Matrix_DIsqrtSym)),              *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
528              *Lx = REAL(GET_SLOT(lme, Matrix_LxSym)),              *Lx = REAL(GET_SLOT(x, Matrix_LxSym)),
529              *RXX = REAL(GET_SLOT(lme, Matrix_RXXSym)),              *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
530              *RZX = REAL(GET_SLOT(lme, Matrix_RZXSym)),              *RZX = REAL(GET_SLOT(x, Matrix_RZXSym)),
531              *dcmp = REAL(getAttrib(lme, Matrix_devCompSym)),              *dcmp = REAL(getAttrib(x, Matrix_devCompSym)),
532              *deviance = REAL(getAttrib(lme, Matrix_devianceSym)),              *deviance = REAL(getAttrib(x, Matrix_devianceSym)),
533              minus1 = -1.,              minus1 = -1.,
534              one = 1.;              one = 1.;
535
536          ssclme_inflate_and_factor(lme);          ssclme_inflate_and_factor(x);
537                                  /* Accumulate logdet of ZtZ+W */                                  /* Accumulate logdet of ZtZ+W */
538          dcmp[0] = dcmp[1] = dcmp[2] = dcmp[3] = 0.;          dcmp[0] = dcmp[1] = dcmp[2] = dcmp[3] = 0.;
539          for (i = 0; i < n; i++) dcmp[0] += log(D[i]);          for (i = 0; i < n; i++) dcmp[0] += log(D[i]);
# Line 579  Line 563
563              }              }
564          }          }
565                                  /* ldl_lsolve on Z'X */                                  /* ldl_lsolve on Z'X */
566          Memcpy(RZX, REAL(GET_SLOT(lme, Matrix_ZtXSym)), n * pp1);          Memcpy(RZX, REAL(GET_SLOT(x, Matrix_ZtXSym)), n * pp1);
567          for (i = 0; i < pp1; i++) {          for (i = 0; i < pp1; i++) {
568              int j;              int j;
569              double *RZXi = RZX + i * n;              double *RZXi = RZX + i * n;
# Line 587  Line 571
571              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
572          }          }
573                                  /* downdate and factor X'X */                                  /* downdate and factor X'X */
574          Memcpy(RXX, REAL(GET_SLOT(lme, Matrix_XtXSym)), pp1 * pp1);          Memcpy(RXX, REAL(GET_SLOT(x, Matrix_XtXSym)), pp1 * pp1);
575          F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,          F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,
576                          RZX, &n, &one, RXX, &pp1);                          RZX, &n, &one, RXX, &pp1);
577          F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);          F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);
578          if (i)          if (i) {
579              error("DPOTRF returned error code %d", i);              warning("Could not factor downdated X'X, code %d", i);
580                dcmp[2] = dcmp[3] = deviance[0] = deviance[1] = NA_REAL;
581            } else {
582                                  /* logdet of RXX */                                  /* logdet of RXX */
583          for (i = 0; i < (pp1 - 1); i++)          for (i = 0; i < (pp1 - 1); i++)
584              dcmp[2] += 2 * log(RXX[i*pp2]);              dcmp[2] += 2 * log(RXX[i*pp2]);
# Line 602  Line 588
588              dcmp[0] - dcmp[1] + nobs*(1+dcmp[3]+log(2*PI/nobs));              dcmp[0] - dcmp[1] + nobs*(1+dcmp[3]+log(2*PI/nobs));
589          deviance[1] = dcmp[0] - dcmp[1] + /* REML */          deviance[1] = dcmp[0] - dcmp[1] + /* REML */
590              dcmp[2] + nreml * (1. + dcmp[3] + log(2. * PI/nreml));              dcmp[2] + nreml * (1. + dcmp[3] + log(2. * PI/nreml));
591            }
592          status[0] = 1;          status[0] = 1;
593          status[1] = 0;          status[1] = 0;
594      }      }
595      return R_NilValue;      return R_NilValue;
596  }  }
597
598    /**
599     * Return the position of probe in the sorted index vector ind.  It is
600     * known that the position is greater than or equal to start so a linear
601     * search from start is used.
602     *
603     * @param probe value to be matched
604     * @param start index at which to start
605     * @param ind vector of indices
606     *
607     * @return index of the entry matching probe
608     */
609  static  static
610  int ldl_update_ind(int probe, int start, const int ind[])  int ldl_update_ind(int probe, int start, const int ind[])
611  {  {
# Line 617  Line 615
615  }  }
616
617  /**  /**
618   * Create the inverse of L and update the diagonal blocks of the inverse   * Update the diagonal blocks of the inverse of LDL' (=Z'Z+W).  The
619   * of LDL' (=Z'Z+W)   * lower Cholesky factors of the updated blocks are stored in the bVar
620     * slot.
621   *   *
622   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
623   *   *
624   * @return R_NilValue (x is updated in place)   * @return R_NilValue (x is updated in place)
625
626   */   */
627    static
628  SEXP ldl_inverse(SEXP x)  SEXP ldl_inverse(SEXP x)
629  {  {
630      SEXP      SEXP
631          Gpsl = GET_SLOT(x, Matrix_GpSym),          Gpsl = GET_SLOT(x, Matrix_GpSym),
LIisl = GET_SLOT(x, Matrix_LIiSym),
LIpsl = GET_SLOT(x, Matrix_LIpSym),
632          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
633      int *Gp = INTEGER(Gpsl),      int *Gp = INTEGER(Gpsl),
634          *Li,          *Parent = INTEGER(GET_SLOT(x, Matrix_ParentSym)),
635          *LIp = INTEGER(LIpsl), *Lp,          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
636          i,          i,
637          nf = length(Gpsl) - 1,          nf = length(Gpsl) - 1,
638          nzc = length(LIpsl) - 1;          nzc = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];
639      double      int maxod = Parent[nzc];
640          *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),      double *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym));
*Lx;
641
642      ssclme_factor(x);      ssclme_factor(x);
643      if (LIp[nzc] == 0) {        /* L and LI are the identity */      if (maxod == 0) {           /* L and L^{-1} are the identity */
644          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
645              Memcpy(REAL(VECTOR_ELT(bVar, i)), DIsqrt + Gp[i],              Memcpy(REAL(VECTOR_ELT(bVar, i)), DIsqrt + Gp[i],
646                     Gp[i+1] - Gp[i]);                     Gp[i+1] - Gp[i]);
647          }          }
648          return R_NilValue;      } else {
649      }          int *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
650      Lp = INTEGER(GET_SLOT(x, Matrix_LpSym));              *Li = INTEGER(GET_SLOT(x, Matrix_LiSym));
Li = INTEGER(GET_SLOT(x, Matrix_LiSym));
Lx = REAL(GET_SLOT(x, Matrix_LxSym));
if (length(LIisl) == LIp[nzc]) { /* LIi is filled */
int *LIi = INTEGER(LIisl),
*nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
j, jj, k, kk, p1, p2, pi1, pi2;

double *LIx = REAL(GET_SLOT(x, Matrix_LIxSym)),
one = 1., zero = 0.;
651
652          memset(LIx, 0, sizeof(double) * LIp[nzc]);          double one = 1.0, zero = 0.,
653                                  /* calculate inverse */              *Lx = REAL(GET_SLOT(x, Matrix_LxSym));
for (i = 0; i < nzc; i++) {
p1 = Lp[i]; p2 = Lp[i+1]; pi1 = LIp[i]; pi2 = LIp[i+1];
/* initialize from unit diagonal term */
kk = pi1;
for (j = p1; j < p2; j++) {
k = Li[j];
while (LIi[kk] < k && kk < pi2) kk++;
if (LIi[kk] != k) error("logic error in ldl_inverse");
LIx[kk] = -Lx[j];
}
for (j = pi1; j < pi2; j++) {
jj = LIi[j];
p1 = Lp[jj]; p2 = Lp[jj+1];
kk = j;
for (jj = p1; jj < p2; jj++) {
k = Li[jj];
while (LIi[kk] < k && kk < pi2) kk++;
if (LIi[kk] != k) error("logic error in ldl_inverse");
LIx[kk] -= Lx[jj]*LIx[j];
}
}
}
for (i = 0; i < nf; i++) { /* accumulate bVar */
int G1 = Gp[i], G2 = Gp[i+1], j, k, kk,
nci = nc[i], nr, nr1, rr;
double *bVi = REAL(VECTOR_ELT(bVar, i)), *tmp;

nr = -1;
for (j = G1; j < G2; j += nci) {
rr = 1 + LIp[j + 1] - LIp[j];
if (rr > nr) nr = rr;
}
tmp = Calloc(nr * nci, double); /* scratch storage */
nr1 = nr + 1;
/* initialize bVi to zero (cosmetic) */
memset(bVi, 0, sizeof(double) * (G2 - G1) * nci);
for (j = G1; j < G2; j += nci) {
memset(tmp, 0, sizeof(double) * nr * nci);
rr = 1 + LIp[j + 1] - LIp[j];
for (k = 0; k < nci; k++) { /* copy columns */
tmp[k * nr1] = 1.; /* (unstored) diagonal elt  */
Memcpy(tmp + k*nr1 + 1, LIx + LIp[j + k], rr - k - 1);
}
/* scale the rows */
tmp[0] = DIsqrt[j]; /* first row only has one non-zero */
for (kk = 1; kk < rr; kk++) {
for (k = 0; k < nci; k++) {
tmp[k * nr + kk] *= DIsqrt[LIi[LIp[j] + kk - 1]];
}
}
F77_CALL(dsyrk)("U", "T", &nci, &rr, &one, tmp, &nr,
&zero, bVi + (j - G1) * nci, &nci);
F77_CALL(dpotrf)("U", &nci, bVi + (j - G1) * nci,
&nci, &kk);
if (kk)         /* should never happen */
error(
"Rank deficient variance matrix at group %d, level %d",
i + 1, j + 1);
}
}
return R_NilValue;
}
if (length(LIisl)) error("logic error in ssclme_ldl_inverse");
else {                  /* LIi and LIx are too big and not used */
int *counts = Calloc(nzc, int), info, maxod = -1;
int *Parent = INTEGER(GET_SLOT(x, Matrix_ParentSym));
int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym));
double one = 1.0, zero = 0.;
/* determine maximum # of off-diagonals */
for (i = nzc - 1; i >= 0; i--) { /* in a column of L^{-1} */
counts[i] = (Parent[i] < 0) ? 0 : 1 + counts[Parent[i]];
if (counts[i] > maxod) maxod = counts[i];
}
Free(counts);
654
655          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
656              int j, jj, k, kk, nci = nc[i], nr, p, p2, pp,              int j, jj, k, kk, nci = nc[i], nr, p, p2, pj, pp,
657                  m = maxod + nci,                  m = maxod + 1,
658                  *ind = Calloc(m, int);                  *ind = Calloc(m, int), G1 = Gp[i], G2 = Gp[i+1];
659              double              double
660                  *tmp = Calloc(m * nci, double),                  *tmp = Calloc(m * nci, double),
661                  *mpt = REAL(VECTOR_ELT(bVar, i));                  *bVi = REAL(VECTOR_ELT(bVar, i));
662
663              for (j = Gp[i]; j < Gp[i+1]; j += nci) {                                  /* initialize bVi to zero */
664                  memset((void *) tmp, 0, sizeof(double) * m * nci);              memset(bVi, 0, sizeof(double) * (G2 - G1) * nci);
665
666                  kk = 0;         /* ind holds indices of non-zeros */              for (j = G1; j < G2; j += nci) {
667                    kk = 0;         /* ind gets indices of non-zeros */
668                  jj = j;         /* in this block of columns */                  jj = j;         /* in this block of columns */
669                  while (jj >= 0) {                  while (jj >= 0) {
670                      ind[kk++] = jj;                      ind[kk++] = jj;
# Line 761  Line 676
676                  for (k = 0; k < nci; k++) {                  for (k = 0; k < nci; k++) {
677                      double *ccol = tmp + k * nr;                      double *ccol = tmp + k * nr;
678
679                      ccol[k] = 1.;                      for (kk = 0; kk < nr; kk++) ccol[kk] = 0.;
680                      kk = k;                      ccol[k] = 1.; /* initialize from unit diagonal */
681                      for (jj = j + k; jj >= 0; jj = Parent[jj]) {                      for (jj = j + k; jj >= 0; jj = Parent[jj]) {
682                          p2 = Lp[jj+1];                          p2 = Lp[jj+1];
683                          pp = kk;                          pp = pj = ldl_update_ind(jj, 0, ind);
684                          for (p = Lp[jj]; p < p2; p++) {                          for (p = Lp[jj]; p < p2; p++) {
685                              pp = ldl_update_ind(Li[p], pp, ind);                              pp = ldl_update_ind(Li[p], pp, ind);
686                              ccol[pp] -= Lx[p] * ccol[kk];                              ccol[pp] -= Lx[p] * ccol[pj];
687                          }                          }
688                      }                      }
689                  }                  }
690                                  /* scale rows */
691                  for (kk = 0; kk < nr; kk++) {                  for (kk = 0; kk < nr; kk++) { /* scale rows */
692                      for (k = 0; k < nci; k++) {                      for (k = 0; k < nci; k++) {
693                          tmp[k * nr + kk] *= DIsqrt[ind[kk]];                          tmp[k * nr + kk] *= DIsqrt[ind[kk]];
694                      }                      }
695                  }                  }
696                  F77_CALL(dsyrk)("U", "T", &nci, &nr, &one, tmp, &nr,                  F77_CALL(dsyrk)("L", "T", &nci, &nr, &one, tmp, &nr,
697                                  &zero, mpt + (j - Gp[i])*nci, &nci);                                  &zero, bVi + (j - G1)*nci, &nci);
698                  F77_CALL(dpotrf)("U", &nci, mpt + (j - Gp[i])*nci,                  F77_CALL(dpotrf)("L", &nci, bVi + (j - G1)*nci,
699                                   &nci, &info);                                   &nci, &jj);
700                  if (info)       /* should never happen */                  if (jj)         /* should never happen */
701                      error(                      error(
702                          "Rank deficient variance matrix at group %d, level %d",                          "Rank deficient variance matrix at group %d, level %d, error code %d",
703                          i + 1, j + 1);                          i + 1, j + 1, jj);
704              }              }
705              Free(tmp); Free(ind);              Free(tmp); Free(ind);
706          }          }
# Line 793  Line 708
708      return R_NilValue;      return R_NilValue;
709  }  }
710
711  SEXP ssclme_invert(SEXP lme)  /**
712     * If necessary, factor Z'Z+Omega, ZtX, and XtX then, if necessary,
713     * form RZX, RXX, and bVar for the inverse of the Cholesky factor.
714     *
715     * @param x pointer to an ssclme object
716     *
717     * @return NULL (x is updated in place)
718     */
719    SEXP ssclme_invert(SEXP x)
720  {  {
721      int *status = LOGICAL(GET_SLOT(lme, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
722      if (!status[0]) ssclme_factor(lme);      if (!status[0]) ssclme_factor(x);
723        if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0]))
724            error("Unable to invert singular factor of downdated X'X");
725      if (!status[1]) {      if (!status[1]) {
726          SEXP          SEXP
727              RZXsl = GET_SLOT(lme, Matrix_RZXSym);              RZXsl = GET_SLOT(x, Matrix_RZXSym);
728          int          int
729              *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),              *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
730              *Li = INTEGER(GET_SLOT(lme, Matrix_LiSym)),              *Li = INTEGER(GET_SLOT(x, Matrix_LiSym)),
731              *Lp = INTEGER(GET_SLOT(lme, Matrix_LpSym)),              *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
732              i,              i,
733              n = dims[0],              n = dims[0],
734              pp1 = dims[1];              pp1 = dims[1];
735          double          double
736              *DIsqrt = REAL(GET_SLOT(lme, Matrix_DIsqrtSym)),              *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
737              *Lx = REAL(GET_SLOT(lme, Matrix_LxSym)),              *Lx = REAL(GET_SLOT(x, Matrix_LxSym)),
738              *RXX = REAL(GET_SLOT(lme, Matrix_RXXSym)),              *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
739              *RZX = REAL(RZXsl),              *RZX = REAL(RZXsl),
740              one = 1.;              one = 1.;
741
# Line 824  Line 749
749              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
750              ldl_ltsolve(n, RZXi, Lp, Li, Lx);              ldl_ltsolve(n, RZXi, Lp, Li, Lx);
751          }          }
752          ldl_inverse(lme);          ldl_inverse(x);
753          status[1] = 1;          status[1] = 1;
754      }      }
755      return R_NilValue;      return R_NilValue;
756  }  }
757
758    /**
759     * Create and insert initial values for Omega_i.
760     *
761     * @param x pointer to an ssclme object
762     *
763     * @return NULL
764     */
765  SEXP ssclme_initial(SEXP x)  SEXP ssclme_initial(SEXP x)
766  {  {
767      SEXP Gpsl = GET_SLOT(x, Matrix_GpSym),      SEXP Gpsl = GET_SLOT(x, Matrix_GpSym),
# Line 869  Line 801
801
802  /**  /**
803   * Extract the conditional estimates of the fixed effects   * Extract the conditional estimates of the fixed effects
* FIXME: Add names
804   *   *
805   * @param x Pointer to an ssclme object   * @param x Pointer to an ssclme object
806   *   *
# Line 897  Line 828
828
829  /**  /**
830   * Extract the conditional modes of the random effects.   * Extract the conditional modes of the random effects.
* FIXME: Change the returned value to be a named list of matrices
*        with dimnames.
831   *   *
832   * @param x Pointer to an ssclme object   * @param x Pointer to an ssclme object
833   *   *
# Line 924  Line 853
853      ssclme_invert(x);      ssclme_invert(x);
854      ryyinv = REAL(GET_SLOT(x, Matrix_RXXSym))[pp1*pp1 - 1];      ryyinv = REAL(GET_SLOT(x, Matrix_RXXSym))[pp1*pp1 - 1];
855      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
856          int nci = nc[i], Mi = (Gp[i+1] - Gp[i]), mi = Mi/nci;          int nci = nc[i], Mi = Gp[i+1] - Gp[i];
857          double *mm;          double *mm;
858
859          SET_VECTOR_ELT(val, i, allocMatrix(REALSXP, nci, mi));          SET_VECTOR_ELT(val, i, allocMatrix(REALSXP, nci, Mi/nci));
860          mm = Memcpy(REAL(VECTOR_ELT(val, i)), b, Mi);          mm = Memcpy(REAL(VECTOR_ELT(val, i)), b, Mi);
861          b += nci * mi;          b += Mi;
862          for (j = 0; j < Mi; j++) mm[j] /= ryyinv;          for (j = 0; j < Mi; j++) mm[j] /= ryyinv;
863      }      }
864      UNPROTECT(1);      UNPROTECT(1);
865      return val;      return val;
866  }  }
867
868  /**  /**
869   * Extract the ML or REML conditional estimate of sigma   * Extract the ML or REML conditional estimate of sigma
870   *   *
# Line 947  Line 877
877  {  {
878      SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);      SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);
879      int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1],      int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1],
880          nobs = INTEGER(GET_SLOT(x, Matrix_ncSym))          nobs = INTEGER(GET_SLOT(x, Matrix_ncSym))[
881          [length(GET_SLOT(x, Matrix_GpSym))];              length(GET_SLOT(x, Matrix_GpSym))];
882
883      ssclme_invert(x);      ssclme_invert(x);
884      return ScalarReal(1./(REAL(RXXsl)[pp1*pp1 - 1] *      return ScalarReal(1./(REAL(RXXsl)[pp1*pp1 - 1] *
# Line 956  Line 886
886                                          nobs + 1 - pp1 : nobs))));                                          nobs + 1 - pp1 : nobs))));
887  }  }
888
889    /**
890     * Calculate the length of the parameter vector, which is called coef
891     * for historical reasons.
892     *
893     * @param nf number of factors
894     * @param nc number of columns in the model matrices for each factor
895     *
896     * @return total length of the coefficient vector
897     */
898  static  static
899  int coef_length(int nf, const int nc[])  int coef_length(int nf, const int nc[])
900  {  {
# Line 965  Line 904
904  }  }
905
906  /**  /**
907   * Extract the upper triangles of the Omega matrices.   * Extract the upper triangles of the Omega matrices.  These aren't
908   * (These are not in any sense "coefficients" but the extractor is   * "coefficients" but the extractor is called coef for historical
909   * called coef for historical reasons.)   * reasons.  Within each group these values are in the order of the
910     * diagonal entries first then the strict upper triangle in row
911     * order.
912   *   *
913   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
914   *   *
# Line 984  Line 925
925
926      vind = 0;      vind = 0;
927      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
928          int j, k, nci = nc[i];          int nci = nc[i];
929            if (nci == 1) {
930                vv[vind++] = REAL(VECTOR_ELT(Omega, i))[0];
931            } else {
932                int j, k, odind = vind + nci, ncip1 = nci + 1;
933          double *omgi = REAL(VECTOR_ELT(Omega, i));          double *omgi = REAL(VECTOR_ELT(Omega, i));
934
935          for (j = 0; j < nci; j++) {          for (j = 0; j < nci; j++) {
936              for (k = 0; k <= j; k++) {                  vv[vind++] = omgi[j * ncip1];
937                  vv[vind++] = omgi[j*nci + k];                  for (k = j + 1; k < nci; k++) {
938                        vv[odind++] = omgi[k*nci + j];
939                    }
940              }              }
941                vind = odind;
942          }          }
943      }      }
944      UNPROTECT(1);      UNPROTECT(1);
# Line 997  Line 946
946  }  }
947
948  /**  /**
949   * Extract the upper triangles of the Omega matrices in the unconstrained   * Extract the unconstrained parameters that determine the
950   * parameterization.   * Omega matrices. (Called coef for historical reasons.)  The
951   * (These are not in any sense "coefficients" but the extractor is   * unconstrained parameters are derived from the LDL' decomposition of
952   * called coef for historical reasons.)   * Omega_i.  The first nc[i] entries in each group are the diagonals
953     * of log(D) followed by the strict lower triangle of L in column
954     * order.
955   *   *
956   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
957   *   *
958   * @return numeric vector of the values in the upper triangles of the   * @return numeric vector of unconstrained parameters that determine the
959   * Omega matrices   * Omega matrices
960   */   */
961  SEXP ssclme_coefUnc(SEXP x)  SEXP ssclme_coefUnc(SEXP x)
# Line 1026  Line 977
977                                   REAL(VECTOR_ELT(Omega, i)), ncisq);                                   REAL(VECTOR_ELT(Omega, i)), ncisq);
978              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
979              if (j)              /* should never happen */              if (j)              /* should never happen */
980                  error("DPOTRF returned error code %d", j);                  error("DPOTRF returned error code %d on Omega[[%d]]",
981                          j, i+1);
982              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
983                  double diagj = tmp[j * ncip1];                  double diagj = tmp[j * ncip1];
984                  vv[vind++] = 2. * log(diagj);                  vv[vind++] = 2. * log(diagj);
# Line 1047  Line 999
999  }  }
1000
1001  /**  /**
1002   * Assign the upper triangles of the Omega matrices in the unconstrained   * Assign the Omega matrices from the unconstrained parameterization.
* parameterization.
1003   *   *
1004   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1005   * @param coef pointer to an numeric vector of appropriate length   * @param coef pointer to an numeric vector of appropriate length
# Line 1079  Line 1030
1030              double              double
1031                  *omgi = REAL(VECTOR_ELT(Omega, i)),                  *omgi = REAL(VECTOR_ELT(Omega, i)),
1032                  *tmp = Calloc(ncisq, double),                  *tmp = Calloc(ncisq, double),
1033                  diagj, one = 1.;                  diagj, one = 1., zero = 0.;
1034                                  /* LD in omgi and L' in tmp */
1035              memset(omgi, 0, sizeof(double) * ncisq);              memset(omgi, 0, sizeof(double) * ncisq);
1036              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
1037                  omgi[j * ncip1] = diagj = exp(cc[cind++]);                  tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);
1038                  for (k = j + 1; k < nci; k++) {                  for (k = j + 1; k < nci; k++) {
1039                      omgi[j*nci + k] = diagj * (tmp[k*nci + j] = cc[odind++]);                      tmp[k*nci + j] = cc[odind++] * diagj;
1040                  }                  }
1041              }              }
1042              F77_CALL(dtrmm)("R", "U", "N", "U", &nci, &nci, &one,              F77_CALL(dsyrk)("U", "T", &nci, &nci, &one,
1043                              tmp, &nci, omgi, &nci);                              tmp, &nci, &zero, omgi, &nci);
1044              Free(tmp);              Free(tmp);
1045              cind = odind;              cind = odind;
1046          }          }
# Line 1100  Line 1051
1051
1052  /**  /**
1053   * Assign the upper triangles of the Omega matrices.   * Assign the upper triangles of the Omega matrices.
1054   * (These are not in any sense "coefficients" but are   * (Called coef for historical reasons.)
* called coef for historical reasons.)
1055   *   *
1056   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1057   * @param coef pointer to an numeric vector of appropriate length   * @param coef pointer to an numeric vector of appropriate length
# Line 1121  Line 1071
1071                coef_length(nf, nc));                coef_length(nf, nc));
1072      cind = 0;      cind = 0;
1073      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1074          int j, k, nci = nc[i];          int nci = nc[i];
1075            if (nci == 1) {
1076                REAL(VECTOR_ELT(Omega, i))[0] = cc[cind++];
1077            } else {
1078                int j, k, odind = cind + nci, ncip1 = nci + 1;
1079          double *omgi = REAL(VECTOR_ELT(Omega, i));          double *omgi = REAL(VECTOR_ELT(Omega, i));
1080
1081          for (j = 0; j < nci; j++) {          for (j = 0; j < nci; j++) {
1082              for (k = 0; k <= j; k++) {                  omgi[j * ncip1] = cc[cind++];
1083                  omgi[j*nci + k] = cc[cind++];                  for (k = j + 1; k < nci; k++) {
1084                        omgi[k*nci + j] = cc[odind++];
1085              }              }
1086          }          }
1087                cind = odind;
1088            }
1089      }      }
1090      status[0] = status[1] = 0;      status[0] = status[1] = 0;
1091      return x;      return x;
1092  }  }
1093
1094    /**
1095     * Perform a number of ECME steps for the REML or ML criterion.
1096     *
1097     * @param x pointer to an ssclme object
1098     * @param nsteps pointer to an integer scalar giving the number of ECME steps to perform
1099     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1100     * @param verb pointer to a logical scalar indicating verbose mode
1101     *
1102     * @return NULL
1103     */
1104  SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)  SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)
1105  {  {
1106      SEXP      SEXP
# Line 1164  Line 1132
1132
1133      p = pp1 - 1;      p = pp1 - 1;
1134      b = RZX + p * n;      b = RZX + p * n;
if (verbose) Rprintf("  EM iterations\n");
for (iter = 0; iter <= nEM; iter++) {
ssclme_invert(x);
1135          if (verbose) {          if (verbose) {
1136              SEXP coef = PROTECT(ssclme_coef(x));              SEXP coef = PROTECT(ssclme_coef(x));
1137              int lc = length(coef); double *cc = REAL(coef);              int lc = length(coef); double *cc = REAL(coef);
1138              Rprintf("%3d %.3f", iter, dev[REML ? 1 : 0]);
1139            ssclme_factor(x);
1140            Rprintf("  EM iterations\n");
1141            Rprintf("%3d %.3f", 0, dev[REML ? 1 : 0]);
1142              for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);              for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);
1143              Rprintf("\n");              Rprintf("\n");
1144              UNPROTECT(1);              UNPROTECT(1);
1145          }          }
1146        for (iter = 0; iter < nEM; iter++) {
1147            ssclme_invert(x);
1148          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
1149              int ki = Gp[i+1] - Gp[i],              int ki = Gp[i+1] - Gp[i],
1150                  nci = nc[i],                  nci = nc[i],
# Line 1191  Line 1161
1161                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,
1162                              &one, vali, &nci);                              &one, vali, &nci);
1163              if (REML) {              if (REML) {
1164                  int mp = mi * p;                  int j;
1165                  F77_CALL(dsyrk)("U", "N", &nci, &mp,                  for (j = 0; j < p; j++) {
1166                                  &alpha, RZX + Gp[i], &nci,                      F77_CALL(dsyrk)("U", "N", &nci, &mi,
1167                                    &alpha, RZX + Gp[i] + j*n, &nci,
1168                                  &one, vali, &nci);                                  &one, vali, &nci);
1169              }              }
1170                }
1171              F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);              F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);
1172              if (info) error("DPOTRF returned error code %d", info);              if (info)
1173                    error("DPOTRF returned error code %d in Omega[[%d]] update",
1174                          info, i + 1);
1175              F77_CALL(dpotri)("U", &nci, vali, &nci, &info);              F77_CALL(dpotri)("U", &nci, vali, &nci, &info);
1176              if (info) error("DPOTRF returned error code %d", info);              if (info)
1177                    error("DPOTRI returned error code %d in Omega[[%d]] update",
1178                          info, i + 1);
1179          }          }
1180          status[0] = status[1] = 0;          status[0] = status[1] = 0;
1181            if (verbose) {
1182                SEXP coef = PROTECT(ssclme_coef(x));
1183                int lc = length(coef); double *cc = REAL(coef);
1184
1185                ssclme_factor(x);
1186                Rprintf("%3d %.3f", iter + 1, dev[REML ? 1 : 0]);
1187                for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);
1188                Rprintf("\n");
1189                UNPROTECT(1);
1190            }
1191      }      }
1192      ssclme_factor(x);      ssclme_factor(x);
1193      return R_NilValue;      return R_NilValue;
1194  }  }
1195
1196  SEXP ssclme_asSscMatrix(SEXP x)  /**
1197     * Return the gradient of the ML or REML deviance.
1198     *
1199     * @param x pointer to an ssclme object
1200     * @param REMLp pointer to a logical scalar indicating if REML is to be used
1201     * @param Uncp pointer to a logical scalar indicating if the unconstrained parameterization is to be used
1202     *
1203     * @return pointer to a numeric vector of the gradient.
1204     */
1205    SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)
1206  {  {
1207      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("sscMatrix")));      SEXP
1208      int *dims = INTEGER(GET_SLOT(val, Matrix_DimSym));          Omega = GET_SLOT(x, Matrix_OmegaSym),
1209            RZXsl = GET_SLOT(x, Matrix_RZXSym),
1210            bVar = GET_SLOT(x, Matrix_bVarSym);
1211        int
1212            *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1213            *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1214            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1215            REML = asLogical(REMLp),
1216            cind, i, n = dims[0],
1217            nf = length(Omega),
1218            nobs, p, pp1 = dims[1],
1219            uncst = asLogical(Uncp);
1220        double
1221            *RZX = REAL(RZXsl),
1222            *b,
1223            alpha,
1224            one = 1.;
1225        SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
1226
1227      dims[0] = dims[1] = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];      ssclme_factor(x);
1228      SET_SLOT(val, Matrix_pSym, duplicate(GET_SLOT(x, Matrix_pSym)));      if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0])) {
1229      SET_SLOT(val, Matrix_iSym, duplicate(GET_SLOT(x, Matrix_iSym)));          int ncoef = coef_length(nf, nc);
1230      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));          for (i = 0; i < ncoef; i++) REAL(ans)[i] = NA_REAL;
CHAR(STRING_ELT(GET_SLOT(x, Matrix_uploSym), 0))[0] = 'U';
1231      UNPROTECT(1);      UNPROTECT(1);
1232      return val;          return ans;
1233  }  }
1234        nobs = nc[nf + 1];
1235        p = pp1 - 1;
1236        b = RZX + p * n;
1237        ssclme_invert(x);
1238        cind = 0;
1239        for (i = 0; i < nf; i++) {
1240            int j, ki = Gp[i+1] - Gp[i],
1241                nci = nc[i], ncip1 = nci + 1, ncisq = nci * nci,
1242                mi = ki/nci;
1243            double
1244                *chol = Memcpy(Calloc(ncisq, double),
1245                               REAL(VECTOR_ELT(Omega, i)), ncisq),
1246                *tmp = Calloc(ncisq, double);
1247
1248  SEXP ssclme_asTscMatrix(SEXP x)
1249            F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1250            if (j)
1251                error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1252            Memcpy(tmp, chol, ncisq);
1253            F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);
1254            if (j)
1255                error("DPOTRI gave error code %d on Omega[[%d]]", j, i + 1);
1256            alpha = (double) -mi;
1257            F77_CALL(dsyrk)("U", "N", &nci, &ki,
1258                            &one, REAL(VECTOR_ELT(bVar, i)), &nci,
1259                            &alpha, tmp, &nci);
1260            alpha = (double)(REML ? (nobs-p) : nobs);
1261            F77_CALL(dsyrk)("U", "N", &nci, &mi,
1262                            &alpha, b + Gp[i], &nci,
1263                            &one, tmp, &nci);
1264            if (REML) {
1265                for (j = 0; j < p; j++) {
1266                    F77_CALL(dsyrk)("U", "N", &nci, &mi,
1267                                    &one, RZX + Gp[i] + j*n, &nci,
1268                                    &one, tmp, &nci);
1269                }
1270            }
1271            if (nci == 1) {
1272                REAL(ans)[cind++] = *tmp *
1273                    (uncst ? *REAL(VECTOR_ELT(Omega, i)) : 1.);
1274            } else {
1275                int k, odind = cind + nci;
1276                if (uncst) {
1277                    int ione = 1, kk;
1278                    double *rr = Calloc(nci, double); /* j'th row of R, the Cholesky factor */
1279                    nlme_symmetrize(tmp, nci);
1280                    for (j = 0; j < nci; j++, cind++) {
1281                        for (k = 0; k < j; k++) rr[k] = 0.;
1282                        for (k = j; k < nci; k++) rr[k] = chol[j + k*nci];
1283                        REAL(ans)[cind] = 0.;
1284                        for (k = j; k < nci; k++) {
1285                            for (kk = j; kk < nci; kk++) {
1286                                REAL(ans)[cind] += rr[k] * rr[kk] *
1287                                    tmp[kk * nci + k];
1288                            }
1289                        }
1290                        for (k = j + 1; k < nci; k++) {
1291                            REAL(ans)[odind++] = 2. * rr[j] *
1292                                F77_CALL(ddot)(&nci, rr, &ione, tmp + k*nci, &ione);
1293                        }
1294                    }
1295                    Free(rr);
1296                } else {
1297                    for (j = 0; j < nci; j++) {
1298                        REAL(ans)[cind++] = tmp[j * ncip1];
1299                        for (k = j + 1; k < nci; k++) {
1300                            REAL(ans)[odind++] = tmp[k*nci + j] * 2.;
1301                        }
1302                    }
1303                }
1304                cind = odind;
1305            }
1306            Free(tmp); Free(chol);
1307        }
1308        UNPROTECT(1);
1309        return ans;
1310    }
1311
1312    /**
1313     * Calculate and return the fitted values.
1314     *
1315     * @param x pointer to an ssclme object
1316     * @param facs list of grouping factors
1317     * @param mmats list of model matrices
1318     * @param useRf pointer to a logical scalar indicating if the random effects should be used
1319     *
1320     * @return pointer to a numeric array of fitted values
1321     */
1322    SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats, SEXP useRf)
1323  {  {
1324      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("tscMatrix")));      SEXP val, b;
1325      int *dims = INTEGER(GET_SLOT(val, Matrix_DimSym));      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1326            i, ione = 1, nf = length(facs), nobs, p;
1327        double *vv, one = 1.0, zero = 0.0;
1328
1329      dims[0] = dims[1] = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];      if (nf < 1)
1330      SET_SLOT(val, Matrix_pSym, duplicate(GET_SLOT(x, Matrix_LpSym)));          error("null factor list passed to ssclme_fitted");
1331      SET_SLOT(val, Matrix_iSym, duplicate(GET_SLOT(x, Matrix_LiSym)));      nobs = length(VECTOR_ELT(facs, 0));
1332      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_LxSym)));      val = PROTECT(allocVector(REALSXP, nobs));
1333      CHAR(STRING_ELT(GET_SLOT(x, Matrix_uploSym), 0))[0] = 'U';      vv = REAL(val);
1334      CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0))[0] = 'U';      p = nc[nf] - 1;
1335        if (p > 0) {
1336            F77_CALL(dgemm)("N", "N", &nobs, &ione, &p, &one,
1337                            REAL(VECTOR_ELT(mmats, nf)), &nobs,
1338                            REAL(PROTECT(ssclme_fixef(x))), &p,
1339                            &zero, vv, &nobs);
1340            UNPROTECT(1);
1341        } else {
1342            memset(vv, 0, sizeof(double) * nobs);
1343        }
1344        if (asLogical(useRf)) {
1345            b = PROTECT(ssclme_ranef(x));
1346            for (i = 0; i < nf; i++) {
1347                int *ff = INTEGER(VECTOR_ELT(facs, i)), j, nci = nc[i];
1348                double *bb = REAL(VECTOR_ELT(b, i)),
1349                    *mm = REAL(VECTOR_ELT(mmats, i));
1350                for (j = 0; j < nobs; ) {
1351                    int nn = 1, lev = ff[j];
1352                    /* check for adjacent rows with same factor level */
1353                    while (ff[j + nn] == lev) nn++;
1354                    F77_CALL(dgemm)("N", "N", &nn, &ione, &nci,
1355                                    &one, mm + j, &nobs,
1356                                    bb + (lev - 1) * nci, &nci,
1357                                    &one, vv + j, &nobs);
1358                    j += nn;
1359                }
1360            }
1361            UNPROTECT(1);
1362        }
1363      UNPROTECT(1);      UNPROTECT(1);
1364      return val;      return val;
1365  }  }
1366
1367    /**
1368     * Return the unscaled variances
1369     *
1370     * @param x pointer to an ssclme object
1371     *
1372     * @return a list similar to the Omega list with the unscaled variances
1373     */
1374    SEXP ssclme_variances(SEXP x)
1375    {
1376        SEXP Omg = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym)));
1377        int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1378            i, nf = length(Omg);
1379
1380        for (i = 0; i < nf; i++) {
1381            double *mm = REAL(VECTOR_ELT(Omg, i));
1382            int j, nci = nc[i];
1383
1384            F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);
1385            if (j)                  /* shouldn't happen */
1386                error("DPOTRF returned error code %d on Omega[%d]",
1387                      j, i + 1);
1388            F77_CALL(dpotri)("U", &nci, mm, &nci, &j);
1389            if (j)                  /* shouldn't happen */
1390                error("DTRTRI returned error code %d on Omega[%d]",
1391                      j, i + 1);
1392            nlme_symmetrize(mm, nci);
1393        }
1394        UNPROTECT(1);
1395        return Omg;
1396    }
1397
1398    /**
1399     * Copy an ssclme object collapsing the fixed effects slots to the response only.
1400     *
1401     * @param x pointer to an ssclme object
1402     *
1403     * @return a duplicate of x with the fixed effects slots collapsed to the response only
1404     */
1405    SEXP ssclme_collapse(SEXP x)
1406    {
1407        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("ssclme"))),
1408            Omega = GET_SLOT(x, Matrix_OmegaSym),
1409            Dim = GET_SLOT(x, Matrix_DimSym);
1410        int nf = length(Omega), nz = INTEGER(Dim)[1];
1411
1412        slot_dup(ans, x, Matrix_DSym);
1413        slot_dup(ans, x, Matrix_DIsqrtSym);
1414        slot_dup(ans, x, Matrix_DimSym);
1415        slot_dup(ans, x, Matrix_GpSym);
1416        slot_dup(ans, x, Matrix_LiSym);
1417        slot_dup(ans, x, Matrix_LpSym);
1418        slot_dup(ans, x, Matrix_LxSym);
1419        slot_dup(ans, x, Matrix_OmegaSym);
1420        slot_dup(ans, x, Matrix_ParentSym);
1421        slot_dup(ans, x, Matrix_bVarSym);
1422        slot_dup(ans, x, Matrix_devianceSym);
1423        slot_dup(ans, x, Matrix_devCompSym);
1424        slot_dup(ans, x, Matrix_iSym);
1425        slot_dup(ans, x, Matrix_ncSym);
1426        slot_dup(ans, x, Matrix_statusSym);
1427        slot_dup(ans, x, Matrix_pSym);
1428        slot_dup(ans, x, Matrix_xSym);
1429        INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;
1430        SET_SLOT(ans, Matrix_XtXSym, allocMatrix(REALSXP, 1, 1));
1431        REAL(GET_SLOT(ans, Matrix_XtXSym))[0] = NA_REAL;
1432        SET_SLOT(ans, Matrix_RXXSym, allocMatrix(REALSXP, 1, 1));
1433        REAL(GET_SLOT(ans, Matrix_RXXSym))[0] = NA_REAL;
1434        SET_SLOT(ans, Matrix_ZtXSym, allocMatrix(REALSXP, nz, 1));
1435        SET_SLOT(ans, Matrix_RZXSym, allocMatrix(REALSXP, nz, 1));
1436        LOGICAL(GET_SLOT(ans, Matrix_statusSym))[0] = 0;
1437        UNPROTECT(1);
1438        return ans;
1439    }
1440
1441
1442    /**
1443     * Create an lme object from its components.  This is not done by
1444     * new("lme", ...) at the R level because of the possibility of
1445     * causing the copying of very large objects.
1446     *
1447     * @param call Pointer to the original call
1448     * @param facs pointer to the list of grouping factors
1449     * @param x pointer to the model matrices (may be of length zero)
1450     * @param model pointer to the model frame
1451     * @param REML pointer to a logical scalar indicating if REML is used
1452     * @param rep pointer to the converged ssclme object
1453     * @param fitted pointer to the fitted values
1454     * @param residuals pointer to the residuals
1455     *
1456     * @return an lme object
1457     */
1458    SEXP ssclme_to_lme(SEXP call, SEXP facs, SEXP x, SEXP model, SEXP REML,
1459                       SEXP rep, SEXP fitted, SEXP residuals)
1460    {
1461        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lme")));
1462
1463        SET_SLOT(ans, install("call"), call);
1464        SET_SLOT(ans, install("facs"), facs);
1465        SET_SLOT(ans, Matrix_xSym, x);
1466        SET_SLOT(ans, install("model"), model);
1467        SET_SLOT(ans, install("REML"), REML);
1468        SET_SLOT(ans, install("rep"), rep);
1469        SET_SLOT(ans, install("fitted"), fitted);
1470        SET_SLOT(ans, install("residuals"), residuals);
1471        UNPROTECT(1);
1472        return ans;
1473    }

Legend:
 Removed from v.22 changed lines Added in v.176

 root@r-forge.r-project.org ViewVC Help Powered by ViewVC 1.0.0
Thanks to: