SCM

SCM Repository

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

Diff of /pkg/src/ssclme.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 100, Sat Apr 17 17:03:59 2004 UTC revision 164, Fri May 14 21:00:04 2004 UTC
# Line 1  Line 1 
1  #include "ssclme.h"  #include "ssclme.h"
2    
 /**  
  * 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;  
 }  
   
 /**  
  * Determine if a fill-reducing permutation is needed for the pairwise  
  * cross-tabulation matrix.  If so, determine such a permutation  
  * (using Metis) then separate the groups.  
  *  
  * @param ctab pointer to a pairwise cross-tabulation object  
  *  
  * @return pointer to an integer R vector.  
  */  
   
 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,  
         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, 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;  
 }  
   
3  static  static
4  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)
5  {  {
# Line 188  Line 99 
99      }      }
100  }  }
101    
102  void ssclme_fill_LIp(int n, const int Parent[], int LIp[])  /**
103     * Calculate and store the maximum number of off-diagonal elements in
104     * the inverse of L, based on the elimination tree.  The maximum is
105     * itself stored in the Parent array.  (FIXME: come up with a better design.)
106     *
107     * @param n number of columns in the matrix
108     * @param Parent elimination tree for the matrix
109     */
110    static void ssclme_calc_maxod(int n, int Parent[])
111  {  {
112      int *sz = Calloc(n, int), i;      int *sz = Calloc(n, int), i, mm = -1;
113      for (i = n - 1; i >= 0; i--) {      for (i = n - 1; i >= 0; i--) {
114          sz[i] = (Parent[i] < 0) ? 0 : 1 + sz[Parent[i]];          sz[i] = (Parent[i] < 0) ? 0 : (1 + sz[Parent[i]]);
115            if (sz[i] > mm) mm = sz[i];
116      }      }
117      LIp[0] = 0;      Parent[n] = mm;
     for (i = 0; i < n; i++) LIp[i+1] = LIp[i] + sz[i];  
118      Free(sz);      Free(sz);
119  }  }
120    
 static  
 void ssclme_fill_LIi(int n, const int Parent[], const int LIp[], int LIi[])  
 {  
     int i;  
     for (i = n; i > 0; i--) {  
         int im1 = i - 1, Par = Parent[im1];  
         if (Par >= 0) {  
             LIi[LIp[im1]] = Par;  
             Memcpy(LIi + LIp[im1] + 1, LIi + LIp[Par],  
                    LIp[Par + 1] - LIp[Par]);  
         }  
     }  
 }  
   
121  SEXP  SEXP
122  ssclme_create(SEXP facs, SEXP ncv, SEXP threshold)  ssclme_create(SEXP facs, SEXP ncv, SEXP threshold)
123  {  {
124      SEXP ctab, nms, ssc, tmp,      SEXP ctab, nms, ssc, tmp,
125          val = PROTECT(allocVector(VECSXP, 2)),          val = PROTECT(allocVector(VECSXP, 2)),
126          dd = PROTECT(allocVector(INTSXP, 3));   /* dimensions of 3-D arrays */          dd = PROTECT(allocVector(INTSXP, 3));   /* dimensions of 3-D arrays */
127      int *Ai, *Ap, *Gp, *LIp, *Lp, *Parent,      int *Ai, *Ap, *Gp, *Lp, *Parent,
128          *nc, Lnz, i, nf = length(facs), nzcol, pp1,          *nc, Lnz, i, nf = length(facs), nzcol, pp1,
129          *dims = INTEGER(dd);          *dims = INTEGER(dd);
130    
# Line 230  Line 135 
135      ssc = VECTOR_ELT(val, 0);      ssc = VECTOR_ELT(val, 0);
136                                  /* Pairwise cross-tabulation */                                  /* Pairwise cross-tabulation */
137      ctab = PROTECT(sscCrosstab(facs, ScalarLogical(1)));      ctab = PROTECT(sscCrosstab(facs, ScalarLogical(1)));
138      SET_VECTOR_ELT(val, 1, ctab_permute(ctab));      SET_VECTOR_ELT(val, 1, sscCrosstab_groupedPerm(ctab));
139      if (length(VECTOR_ELT(val, 1)) > 0) {/* Fill-reducing permutation */      if (length(VECTOR_ELT(val, 1)) > 0) {/* Fill-reducing permutation */
140          ssc_symbolic_permute(INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],          ssc_symbolic_permute(INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],
141                               1, INTEGER(VECTOR_ELT(val, 1)),                               1, INTEGER(VECTOR_ELT(val, 1)),
# Line 252  Line 157 
157      SET_SLOT(ssc, Matrix_RXXSym, allocMatrix(REALSXP, pp1, pp1));      SET_SLOT(ssc, Matrix_RXXSym, allocMatrix(REALSXP, pp1, pp1));
158      SET_SLOT(ssc, Matrix_ZtXSym, allocMatrix(REALSXP, nzcol, pp1));      SET_SLOT(ssc, Matrix_ZtXSym, allocMatrix(REALSXP, nzcol, pp1));
159      SET_SLOT(ssc, Matrix_RZXSym, allocMatrix(REALSXP, nzcol, pp1));      SET_SLOT(ssc, Matrix_RZXSym, allocMatrix(REALSXP, nzcol, pp1));
160          /* Zero the symmetric matrices (for cosmetic reasons only). */                                  /* Zero symmetric matrices (cosmetic) */
161      memset(REAL(GET_SLOT(ssc, Matrix_XtXSym)), 0,      memset(REAL(GET_SLOT(ssc, Matrix_XtXSym)), 0,
162             sizeof(double) * pp1 * pp1);             sizeof(double) * pp1 * pp1);
163      memset(REAL(GET_SLOT(ssc, Matrix_RXXSym)), 0,      memset(REAL(GET_SLOT(ssc, Matrix_RXXSym)), 0,
164             sizeof(double) * pp1 * pp1);             sizeof(double) * pp1 * pp1);
165      SET_SLOT(ssc, Matrix_LpSym, allocVector(INTSXP, nzcol + 1));      SET_SLOT(ssc, Matrix_LpSym, allocVector(INTSXP, nzcol + 1));
166      Lp = INTEGER(GET_SLOT(ssc, Matrix_LpSym));      Lp = INTEGER(GET_SLOT(ssc, Matrix_LpSym));
167      SET_SLOT(ssc, Matrix_ParentSym, allocVector(INTSXP, nzcol));      SET_SLOT(ssc, Matrix_ParentSym, allocVector(INTSXP, nzcol + 1));
168      Parent = INTEGER(GET_SLOT(ssc, Matrix_ParentSym));      Parent = INTEGER(GET_SLOT(ssc, Matrix_ParentSym));
169      SET_SLOT(ssc, Matrix_DSym, allocVector(REALSXP, nzcol));      SET_SLOT(ssc, Matrix_DSym, allocVector(REALSXP, nzcol));
170      SET_SLOT(ssc, Matrix_DIsqrtSym, allocVector(REALSXP, nzcol));      SET_SLOT(ssc, Matrix_DIsqrtSym, allocVector(REALSXP, nzcol));
# Line 267  Line 172 
172                   (int *) R_alloc(nzcol, sizeof(int)), /* Lnz */                   (int *) R_alloc(nzcol, sizeof(int)), /* Lnz */
173                   (int *) R_alloc(nzcol, sizeof(int)), /* Flag */                   (int *) R_alloc(nzcol, sizeof(int)), /* Flag */
174                   (int *) NULL, (int *) NULL); /* P & Pinv */                   (int *) NULL, (int *) NULL); /* P & Pinv */
175        ssclme_calc_maxod(nzcol, Parent);
176      Lnz = Lp[nzcol];      Lnz = Lp[nzcol];
177      SET_SLOT(ssc, Matrix_LiSym, allocVector(INTSXP, Lnz));      SET_SLOT(ssc, Matrix_LiSym, allocVector(INTSXP, Lnz));
178      SET_SLOT(ssc, Matrix_LxSym, allocVector(REALSXP, Lnz));      SET_SLOT(ssc, Matrix_LxSym, allocVector(REALSXP, Lnz));
# Line 304  Line 210 
210          memset(REAL(VECTOR_ELT(tmp, i)), 0,          memset(REAL(VECTOR_ELT(tmp, i)), 0,
211                 sizeof(double) * nci * nci * mi);                 sizeof(double) * nci * nci * mi);
212      }      }
     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);  
     }  
213      UNPROTECT(2);      UNPROTECT(2);
214      return val;      return val;
215  }  }
# Line 369  Line 264 
264          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
265          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
266          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
267            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
268          i, j, k,          i, j, k,
269          ione = 1,          ione = 1,
270          nf = length(mmats) - 1,          nf = length(mmats) - 1,
# Line 432  Line 328 
328          for (k = j+1; k < nf; k++) { /* off-diagonals */          for (k = j+1; k < nf; k++) { /* off-diagonals */
329              int *fpk = INTEGER(VECTOR_ELT(facs, k)),              int *fpk = INTEGER(VECTOR_ELT(facs, k)),
330                  *Apk = Ap + Gp[k],                  *Apk = Ap + Gp[k],
331                  nck = nc[k];                  nck = nc[k],
332                    scalar = ncj == 1 && nck == 1;
333              double              double
334                  *Zk = Z[k];                  *Zk = Z[k], *work;
335                if (!scalar) work = Calloc(ncj * nck, double);
336              for (i = 0; i < nobs; i++) {              for (i = 0; i < nobs; i++) {
337                  int ii, ind = -1, fpji = fpj[i] - 1,                  int ii, ind = -1, fpji = fpj[i] - 1,
338                      row = Gp[j] + fpji * ncj,                      row = Gp[j] + fpji * ncj,
339                      fpki = fpk[i] - 1,                      fpki = fpk[i] - 1,
340                      lastind = Apk[fpki + 1];                      lastind = Apk[fpki*nck + 1];
341                  for (ii = Apk[fpki]; ii < lastind; ii++) {                  for (ii = Apk[fpki*nck]; ii < lastind; ii++) {
342                      if (Ai[ii] == row) {                      if (Ai[ii] == row) {
343                          ind = ii;                          ind = ii;
344                          break;                          break;
345                      }                      }
346                  }                  }
347                  if (ind < 0) error("logic error in ssclme_update_mm");                  if (ind < 0) error("logic error in ssclme_update_mm");
348                  if (Ncj || nck > 1) {                  if (scalar) {   /* update scalars directly */
349                                  /* FIXME: run a loop to update */                      Ax[ind] += Zj[i] * Zk[i];
350                      error("code not yet written");                  } else {
351                  } else {        /* update scalars directly */                      int jj, offset = ind - Apk[fpki * nck];
352                      Ax[ind] += Zj[fpji] * Zk[fpki];                      F77_CALL(dgemm)("T", "N", &ncj, &nck, &ione, &one,
353                                        Zj + i, &nobs, Zk + i, &nobs,
354                                        &zero, work, &ncj);
355                        for (jj = 0; jj < nck; jj++) {
356                            ind = Apk[fpki * nck + jj] + offset;
357                            if (Ai[ind] != row)
358                                error("logic error in ssclme_update_mm");
359                            for (ii = 0; ii < ncj; ii++) {
360                                Ax[ind++] += work[jj * ncj + ii];
361                  }                  }
362              }              }
363          }          }
364      }      }
365                if (!scalar) Free(work);
366            }
367        }
368      Free(Z);      Free(Z);
369      ssclme_transfer_dimnames(x, facs, mmats);      ssclme_transfer_dimnames(x, facs, mmats);
370        status[0] = status[1] = 0;
371      return R_NilValue;      return R_NilValue;
372  }  }
373    
# Line 589  Line 498 
498          F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,          F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,
499                          RZX, &n, &one, RXX, &pp1);                          RZX, &n, &one, RXX, &pp1);
500          F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);          F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);
501          if (i)          if (i) {
502              error("DPOTRF returned error code %d", i);              warning("Could not factor downdated X'X, code %d", i);
503                dcmp[2] = dcmp[3] = deviance[0] = deviance[1] = NA_REAL;
504            } else {
505                                  /* logdet of RXX */                                  /* logdet of RXX */
506          for (i = 0; i < (pp1 - 1); i++)          for (i = 0; i < (pp1 - 1); i++)
507              dcmp[2] += 2 * log(RXX[i*pp2]);              dcmp[2] += 2 * log(RXX[i*pp2]);
# Line 600  Line 511 
511              dcmp[0] - dcmp[1] + nobs*(1+dcmp[3]+log(2*PI/nobs));              dcmp[0] - dcmp[1] + nobs*(1+dcmp[3]+log(2*PI/nobs));
512          deviance[1] = dcmp[0] - dcmp[1] + /* REML */          deviance[1] = dcmp[0] - dcmp[1] + /* REML */
513              dcmp[2] + nreml * (1. + dcmp[3] + log(2. * PI/nreml));              dcmp[2] + nreml * (1. + dcmp[3] + log(2. * PI/nreml));
514            }
515          status[0] = 1;          status[0] = 1;
516          status[1] = 0;          status[1] = 0;
517      }      }
# Line 615  Line 527 
527  }  }
528    
529  /**  /**
530   * 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)
  * of LDL' (=Z'Z+W)  
531   *   *
532   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
533   *   *
# Line 628  Line 539 
539  {  {
540      SEXP      SEXP
541          Gpsl = GET_SLOT(x, Matrix_GpSym),          Gpsl = GET_SLOT(x, Matrix_GpSym),
         LIisl = GET_SLOT(x, Matrix_LIiSym),  
         LIpsl = GET_SLOT(x, Matrix_LIpSym),  
542          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
543      int *Gp = INTEGER(Gpsl),      int *Gp = INTEGER(Gpsl),
544          *Li,          *Parent = INTEGER(GET_SLOT(x, Matrix_ParentSym)),
545          *LIp = INTEGER(LIpsl), *Lp,          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
546          i,          i,
547          nf = length(Gpsl) - 1,          nf = length(Gpsl) - 1,
548          nzc = length(LIpsl) - 1;          nzc = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];
549      double      int maxod = Parent[nzc];
550          *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),      double *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym));
         *Lx;  
551    
552      ssclme_factor(x);      ssclme_factor(x);
553      if (LIp[nzc] == 0) {        /* L and LI are the identity */      if (maxod == 0) {           /* L and L^{-1} are the identity */
554          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
555              Memcpy(REAL(VECTOR_ELT(bVar, i)), DIsqrt + Gp[i],              Memcpy(REAL(VECTOR_ELT(bVar, i)), DIsqrt + Gp[i],
556                     Gp[i+1] - Gp[i]);                     Gp[i+1] - Gp[i]);
557          }          }
558          return R_NilValue;      } else {
559      }          int *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
560      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.;  
   
         memset(LIx, 0, sizeof(double) * LIp[nzc]);  
                                 /* calculate inverse */  
         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;  
561    
562              nr = -1;          double one = 1.0, zero = 0.,
563              for (j = G1; j < G2; j += nci) {              *Lx = REAL(GET_SLOT(x, Matrix_LxSym));
                 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);  
564    
565          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
566              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,
567                  m = maxod + nci,                  m = maxod + 1,
568                  *ind = Calloc(m, int);                  *ind = Calloc(m, int), G1 = Gp[i], G2 = Gp[i+1];
569              double              double
570                  *tmp = Calloc(m * nci, double),                  *tmp = Calloc(m * nci, double),
571                  *mpt = REAL(VECTOR_ELT(bVar, i));                  *bVi = REAL(VECTOR_ELT(bVar, i));
572    
573              for (j = Gp[i]; j < Gp[i+1]; j += nci) {                                  /* initialize bVi to zero */
574                  memset((void *) tmp, 0, sizeof(double) * m * nci);              memset(bVi, 0, sizeof(double) * (G2 - G1) * nci);
575    
576                  kk = 0;         /* ind holds indices of non-zeros */              for (j = G1; j < G2; j += nci) {
577                    kk = 0;         /* ind gets indices of non-zeros */
578                  jj = j;         /* in this block of columns */                  jj = j;         /* in this block of columns */
579                  while (jj >= 0) {                  while (jj >= 0) {
580                      ind[kk++] = jj;                      ind[kk++] = jj;
# Line 760  Line 586 
586                  for (k = 0; k < nci; k++) {                  for (k = 0; k < nci; k++) {
587                      double *ccol = tmp + k * nr;                      double *ccol = tmp + k * nr;
588    
589                      ccol[k] = 1.;                      for (kk = 0; kk < nr; kk++) ccol[kk] = 0.;
590                      kk = k;                      ccol[k] = 1.; /* initialize from unit diagonal */
591                      for (jj = j + k; jj >= 0; jj = Parent[jj]) {                      for (jj = j + k; jj >= 0; jj = Parent[jj]) {
592                          p2 = Lp[jj+1];                          p2 = Lp[jj+1];
593                          pp = kk;                          pp = pj = ldl_update_ind(jj, 0, ind);
594                          for (p = Lp[jj]; p < p2; p++) {                          for (p = Lp[jj]; p < p2; p++) {
595                              pp = ldl_update_ind(Li[p], pp, ind);                              pp = ldl_update_ind(Li[p], pp, ind);
596                              ccol[pp] -= Lx[p] * ccol[kk];                              ccol[pp] -= Lx[p] * ccol[pj];
597                          }                          }
598                      }                      }
599                  }                  }
600                                  /* scale rows */  
601                  for (kk = 0; kk < nr; kk++) {                  for (kk = 0; kk < nr; kk++) { /* scale rows */
602                      for (k = 0; k < nci; k++) {                      for (k = 0; k < nci; k++) {
603                          tmp[k * nr + kk] *= DIsqrt[ind[kk]];                          tmp[k * nr + kk] *= DIsqrt[ind[kk]];
604                      }                      }
605                  }                  }
606                  F77_CALL(dsyrk)("U", "T", &nci, &nr, &one, tmp, &nr,                  F77_CALL(dsyrk)("L", "T", &nci, &nr, &one, tmp, &nr,
607                                  &zero, mpt + (j - Gp[i])*nci, &nci);                                  &zero, bVi + (j - G1)*nci, &nci);
608                  F77_CALL(dpotrf)("U", &nci, mpt + (j - Gp[i])*nci,                  F77_CALL(dpotrf)("L", &nci, bVi + (j - G1)*nci,
609                                   &nci, &info);                                   &nci, &jj);
610                  if (info)       /* should never happen */                  if (jj)         /* should never happen */
611                      error(                      error(
612                          "Rank deficient variance matrix at group %d, level %d",                          "Rank deficient variance matrix at group %d, level %d, error code %d",
613                          i + 1, j + 1);                          i + 1, j + 1, jj);
614              }              }
615              Free(tmp); Free(ind);              Free(tmp); Free(ind);
616          }          }
617      }      }
618      return R_NilValue;      return R_NilValue;
619  }  }
620    
621  SEXP ssclme_invert(SEXP x)  SEXP ssclme_invert(SEXP x)
622  {  {
623      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
624      if (!status[0]) ssclme_factor(x);      if (!status[0]) ssclme_factor(x);
625        if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0]))
626            error("Unable to invert singular factor of downdated X'X");
627      if (!status[1]) {      if (!status[1]) {
628          SEXP          SEXP
629              RZXsl = GET_SLOT(x, Matrix_RZXSym);              RZXsl = GET_SLOT(x, Matrix_RZXSym);
# Line 867  Line 696 
696    
697  /**  /**
698   * Extract the conditional estimates of the fixed effects   * Extract the conditional estimates of the fixed effects
  * FIXME: Add names  
699   *   *
700   * @param x Pointer to an ssclme object   * @param x Pointer to an ssclme object
701   *   *
# Line 895  Line 723 
723    
724  /**  /**
725   * 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.  
726   *   *
727   * @param x Pointer to an ssclme object   * @param x Pointer to an ssclme object
728   *   *
# Line 946  Line 772 
772  {  {
773      SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);      SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);
774      int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1],      int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1],
775          nobs = INTEGER(GET_SLOT(x, Matrix_ncSym))          nobs = INTEGER(GET_SLOT(x, Matrix_ncSym))[
776          [length(GET_SLOT(x, Matrix_GpSym))];              length(GET_SLOT(x, Matrix_GpSym))];
777    
778      ssclme_invert(x);      ssclme_invert(x);
779      return ScalarReal(1./(REAL(RXXsl)[pp1*pp1 - 1] *      return ScalarReal(1./(REAL(RXXsl)[pp1*pp1 - 1] *
# Line 965  Line 791 
791    
792  /**  /**
793   * Extract the upper triangles of the Omega matrices.   * Extract the upper triangles of the Omega matrices.
794   * (These are not in any sense "coefficients" but the extractor is   * (These aren't "coefficients" but the extractor is
795   * called coef for historical reasons.)   * called coef for historical reasons.)
796   *   *
797   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
# Line 983  Line 809 
809    
810      vind = 0;      vind = 0;
811      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
812          int j, k, nci = nc[i];          int nci = nc[i];
813            if (nci == 1) {
814                vv[vind++] = REAL(VECTOR_ELT(Omega, i))[0];
815            } else {
816                int j, k, odind = vind + nci, ncip1 = nci + 1;
817          double *omgi = REAL(VECTOR_ELT(Omega, i));          double *omgi = REAL(VECTOR_ELT(Omega, i));
818    
819          for (j = 0; j < nci; j++) {          for (j = 0; j < nci; j++) {
820              for (k = 0; k <= j; k++) {                  vv[vind++] = omgi[j * ncip1];
821                  vv[vind++] = omgi[j*nci + k];                  for (k = j + 1; k < nci; k++) {
822                        vv[odind++] = omgi[k*nci + j];
823                    }
824              }              }
825                vind = odind;
826          }          }
827      }      }
828      UNPROTECT(1);      UNPROTECT(1);
# Line 996  Line 830 
830  }  }
831    
832  /**  /**
833   * Extract the upper triangles of the Omega matrices in the unconstrained   * Extract the unconstrained parameters that determine the
834   * parameterization.   * Omega matrices. (Called coef for historical reasons.)
  * (These are not in any sense "coefficients" but the extractor is  
  * called coef for historical reasons.)  
835   *   *
836   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
837   *   *
838   * @return numeric vector of the values in the upper triangles of the   * @return numeric vector of unconstrained parameters that determine the
839   * Omega matrices   * Omega matrices
840   */   */
841  SEXP ssclme_coefUnc(SEXP x)  SEXP ssclme_coefUnc(SEXP x)
# Line 1025  Line 857 
857                                   REAL(VECTOR_ELT(Omega, i)), ncisq);                                   REAL(VECTOR_ELT(Omega, i)), ncisq);
858              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
859              if (j)              /* should never happen */              if (j)              /* should never happen */
860                  error("DPOTRF returned error code %d", j);                  error("DPOTRF returned error code %d on Omega[[%d]]",
861                          j, i+1);
862              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
863                  double diagj = tmp[j * ncip1];                  double diagj = tmp[j * ncip1];
864                  vv[vind++] = 2. * log(diagj);                  vv[vind++] = 2. * log(diagj);
# Line 1046  Line 879 
879  }  }
880    
881  /**  /**
882   * Assign the upper triangles of the Omega matrices in the unconstrained   * Assign the Omega matrices from the unconstrained parameterization.
  * parameterization.  
883   *   *
884   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
885   * @param coef pointer to an numeric vector of appropriate length   * @param coef pointer to an numeric vector of appropriate length
# Line 1078  Line 910 
910              double              double
911                  *omgi = REAL(VECTOR_ELT(Omega, i)),                  *omgi = REAL(VECTOR_ELT(Omega, i)),
912                  *tmp = Calloc(ncisq, double),                  *tmp = Calloc(ncisq, double),
913                  diagj, one = 1.;                  diagj, one = 1., zero = 0.;
914              /* FIXEME: Change this to use a factor and dsyrk */  
                                 /* LD in omgi and L' in tmp */  
915              memset(omgi, 0, sizeof(double) * ncisq);              memset(omgi, 0, sizeof(double) * ncisq);
916              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
917                  omgi[j * ncip1] = diagj = exp(cc[cind++]);                  tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);
918                  for (k = j + 1; k < nci; k++) {                  for (k = j + 1; k < nci; k++) {
919                      omgi[j*nci + k] = diagj * (tmp[k*nci + j] = cc[odind++]);                      tmp[k*nci + j] = cc[odind++] * diagj;
920                  }                  }
921              }              }
922              F77_CALL(dtrmm)("R", "U", "N", "U", &nci, &nci, &one,              F77_CALL(dsyrk)("U", "T", &nci, &nci, &one,
923                              tmp, &nci, omgi, &nci);                              tmp, &nci, &zero, omgi, &nci);
924              Free(tmp);              Free(tmp);
925              cind = odind;              cind = odind;
926          }          }
# Line 1100  Line 931 
931    
932  /**  /**
933   * Assign the upper triangles of the Omega matrices.   * Assign the upper triangles of the Omega matrices.
934   * (These are not in any sense "coefficients" but are   * (Called coef for historical reasons.)
  * called coef for historical reasons.)  
935   *   *
936   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
937   * @param coef pointer to an numeric vector of appropriate length   * @param coef pointer to an numeric vector of appropriate length
# Line 1121  Line 951 
951                coef_length(nf, nc));                coef_length(nf, nc));
952      cind = 0;      cind = 0;
953      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
954          int j, k, nci = nc[i];          int nci = nc[i];
955            if (nci == 1) {
956                REAL(VECTOR_ELT(Omega, i))[0] = cc[cind++];
957            } else {
958                int j, k, odind = cind + nci, ncip1 = nci + 1;
959          double *omgi = REAL(VECTOR_ELT(Omega, i));          double *omgi = REAL(VECTOR_ELT(Omega, i));
960    
961          for (j = 0; j < nci; j++) {          for (j = 0; j < nci; j++) {
962              for (k = 0; k <= j; k++) {                  omgi[j * ncip1] = cc[cind++];
963                  omgi[j*nci + k] = cc[cind++];                  for (k = j + 1; k < nci; k++) {
964                        omgi[k*nci + j] = cc[odind++];
965                    }
966              }              }
967                cind = odind;
968          }          }
969      }      }
970      status[0] = status[1] = 0;      status[0] = status[1] = 0;
# Line 1199  Line 1037 
1037                  }                  }
1038              }              }
1039              F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);              F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);
1040              if (info) error("DPOTRF returned error code %d", info);              if (info)
1041                    error("DPOTRF returned error code %d in Omega[[%d]] update",
1042                          info, i + 1);
1043              F77_CALL(dpotri)("U", &nci, vali, &nci, &info);              F77_CALL(dpotri)("U", &nci, vali, &nci, &info);
1044              if (info) error("DPOTRI returned error code %d", info);              if (info)
1045                    error("DPOTRI returned error code %d in Omega[[%d]] update",
1046                          info, i + 1);
1047          }          }
1048          status[0] = status[1] = 0;          status[0] = status[1] = 0;
1049      }      }
# Line 1209  Line 1051 
1051      return R_NilValue;      return R_NilValue;
1052  }  }
1053    
1054  SEXP ssclme_gradient(SEXP x, SEXP REMLp)  SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)
1055  {  {
1056      SEXP      SEXP
1057            Omega = GET_SLOT(x, Matrix_OmegaSym),
1058          RZXsl = GET_SLOT(x, Matrix_RZXSym),          RZXsl = GET_SLOT(x, Matrix_RZXSym),
         ans = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym))),  
1059          ncsl = GET_SLOT(x, Matrix_ncSym),          ncsl = GET_SLOT(x, Matrix_ncSym),
1060          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
1061      int      int
# Line 1221  Line 1063 
1063          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1064          *nc = INTEGER(ncsl),          *nc = INTEGER(ncsl),
1065          REML = asLogical(REMLp),          REML = asLogical(REMLp),
1066          i, info,          cind, i, n = dims[0],
1067          n = dims[0],          nf = length(Omega),
1068          nf = length(ncsl) - 2,          nobs, p, pp1 = dims[1],
1069          nobs = nc[nf + 1],          uncst = asLogical(Uncp);
         p,  
         pp1 = dims[1];  
1070      double      double
1071          *RZX = REAL(RZXsl),          *RZX = REAL(RZXsl),
1072          *b,          *b,
1073          alpha,          alpha,
1074          one = 1.;          one = 1.;
1075        SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
1076    
1077        ssclme_factor(x);
1078        if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0])) {
1079            int ncoef = coef_length(nf, nc);
1080            for (i = 0; i < ncoef; i++) REAL(ans)[i] = NA_REAL;
1081            UNPROTECT(1);
1082            return ans;
1083        }
1084        nobs = nc[nf + 1];
1085      p = pp1 - 1;      p = pp1 - 1;
1086      b = RZX + p * n;      b = RZX + p * n;
1087      ssclme_invert(x);      ssclme_invert(x);
1088        cind = 0;
1089      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1090          int ki = Gp[i+1] - Gp[i],          int j, ki = Gp[i+1] - Gp[i],
1091              nci = nc[i],              nci = nc[i], ncip1 = nci + 1, ncisq = nci * nci,
1092              mi = ki/nci;              mi = ki/nci;
1093          double          double
1094              *vali = REAL(VECTOR_ELT(ans, i));              *chol = Memcpy(Calloc(ncisq, double),
1095                               REAL(VECTOR_ELT(Omega, i)), ncisq),
1096                *tmp = Calloc(ncisq, double);
1097    
1098          F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);  
1099          if (info)          F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1100              error("DPOTRF returned error code %d for component %d of Omega",          if (j)
1101                    info, i + 1);              error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1102          F77_CALL(dpotri)("U", &nci, vali, &nci, &info);          Memcpy(tmp, chol, ncisq);
1103          if (info)          F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);
1104              error("DPOTRI returned error code %d for component %d of Omega",          if (j)
1105                    info, i + 1);              error("DPOTRI gave error code %d on Omega[[%d]]", j, i + 1);
1106          alpha = (double) -mi;          alpha = (double) -mi;
1107          F77_CALL(dsyrk)("U", "N", &nci, &ki,          F77_CALL(dsyrk)("U", "N", &nci, &ki,
1108                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,
1109                          &alpha, vali, &nci);                          &alpha, tmp, &nci);
1110          alpha = ((double)(REML?(nobs-p):nobs));          alpha = ((double)(REML?(nobs-p):nobs));
1111          F77_CALL(dsyrk)("U", "N", &nci, &mi,          F77_CALL(dsyrk)("U", "N", &nci, &mi,
1112                          &alpha, b + Gp[i], &nci,                          &alpha, b + Gp[i], &nci,
1113                          &one, vali, &nci);                          &one, tmp, &nci);
1114          if (REML) {          if (REML) {
             int j;  
1115              for (j = 0; j < p; j++) {              for (j = 0; j < p; j++) {
1116                  F77_CALL(dsyrk)("U", "N", &nci, &mi,                  F77_CALL(dsyrk)("U", "N", &nci, &mi,
1117                                  &one, RZX + Gp[i] + j*n, &nci,                                  &one, RZX + Gp[i] + j*n, &nci,
1118                                  &one, vali, &nci);                                  &one, tmp, &nci);
1119                }
1120            }
1121            if (nci == 1) {
1122                REAL(ans)[cind++] = *tmp *
1123                    (uncst ? *REAL(VECTOR_ELT(Omega, i)) : 1.);
1124            } else {
1125                int k, odind = cind + nci;
1126                if (uncst) {
1127                    int ione = 1, kk;
1128                    double *rr = Calloc(nci, double);
1129                    nlme_symmetrize(tmp, nci);
1130                    for (j = 0; j < nci; j++, cind++) {
1131                        for (k = 0; k < nci; k++) rr[k] = chol[j + k*nci];
1132                        REAL(ans)[cind] = 0.;
1133                        for (k = j; k < nci; k++) {
1134                            for (kk = j; kk < nci; kk++) {
1135                                REAL(ans)[cind] += rr[k] * rr[kk] *
1136                                    tmp[kk * nci + k];
1137                            }
1138                        }
1139                        for (k = 0; k < nci; k++) rr[k] *= rr[j];
1140                        for (k = j + 1; k < nci; k++) {
1141                            REAL(ans)[odind++] =
1142                                F77_CALL(ddot)(&nci, rr, &ione, tmp + k, &nci) +
1143                                F77_CALL(ddot)(&nci, rr, &ione,
1144                                               tmp + k*nci, &ione);
1145                        }
1146                    }
1147                    Free(rr);
1148                } else {
1149                    for (j = 0; j < nci; j++) {
1150                        REAL(ans)[cind++] = tmp[j * ncip1];
1151                        for (k = j + 1; k < nci; k++) {
1152                            REAL(ans)[odind++] = tmp[k*nci + j] * 2.;
1153              }              }
1154          }          }
1155      }      }
1156                cind = odind;
1157            }
1158            Free(tmp); Free(chol);
1159        }
1160      UNPROTECT(1);      UNPROTECT(1);
1161      return ans;      return ans;
1162  }  }
1163    
1164  SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats)  SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats, SEXP useRf)
1165  {  {
1166      SEXP val, b;      SEXP val, b;
1167      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
# Line 1294  Line 1183 
1183      } else {      } else {
1184          memset(vv, 0, sizeof(double) * nobs);          memset(vv, 0, sizeof(double) * nobs);
1185      }      }
1186        if (asLogical(useRf)) {
1187      b = PROTECT(ssclme_ranef(x));      b = PROTECT(ssclme_ranef(x));
1188      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1189          int *ff = INTEGER(VECTOR_ELT(facs, i)), j, nci = nc[i];          int *ff = INTEGER(VECTOR_ELT(facs, i)), j, nci = nc[i];
# Line 1310  Line 1200 
1200              j += nn;              j += nn;
1201          }          }
1202      }      }
1203      UNPROTECT(2);          UNPROTECT(1);
1204        }
1205        UNPROTECT(1);
1206      return val;      return val;
1207  }  }
1208    
1209  SEXP ssclme_variances(SEXP x, SEXP REML)  /**
1210     * Return the unscaled variances
1211     *
1212     * @param x pointer to an ssclme object
1213     *
1214     * @return a list similar to the Omega list with the unscaled variances
1215     */
1216    SEXP ssclme_variances(SEXP x)
1217  {  {
1218      SEXP Omg = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym))),      SEXP Omg = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym)));
         val = PROTECT(allocVector(VECSXP, 2));  
1219      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1220          i, nf = length(Omg);          i, nf = length(Omg);
     double sigmasq;  
   
1221    
     SET_VECTOR_ELT(val, 0, Omg);  
     SET_VECTOR_ELT(val, 1, ssclme_sigma(x, REML));  
     sigmasq = REAL(VECTOR_ELT(val, 1))[0];  
     sigmasq = (sigmasq) * (sigmasq);  
1222      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1223          double *mm = REAL(VECTOR_ELT(Omg, i));          double *mm = REAL(VECTOR_ELT(Omg, i));
1224          int j, k, nci = nc[i], ncip1 = nci+1;          int j, nci = nc[i];
1225    
1226          F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);          F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);
1227          if (j)                  /* shouldn't happen */          if (j)                  /* shouldn't happen */
# Line 1339  Line 1231 
1231          if (j)                  /* shouldn't happen */          if (j)                  /* shouldn't happen */
1232              error("DTRTRI returned error code %d on Omega[%d]",              error("DTRTRI returned error code %d on Omega[%d]",
1233                    j, i + 1);                    j, i + 1);
1234          for (j = 0; j < nci; j++) {          nlme_symmetrize(mm, nci);
             mm[j * ncip1] *= sigmasq;  
             for (k = 0; k < j; k++) {  
                 mm[j + k * nci] = (mm[k + j * nci] *= sigmasq);  
             }  
1235          }          }
1236        UNPROTECT(1);
1237        return Omg;
1238      }      }
1239      UNPROTECT(2);  
1240      return val;  #define slot_dup(sym)  SET_SLOT(ans, sym, duplicate(GET_SLOT(x, sym)))
1241    
1242    SEXP ssclme_collapse(SEXP x)
1243    {
1244        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("ssclme"))),
1245            Omega = GET_SLOT(x, Matrix_OmegaSym),
1246            Dim = GET_SLOT(x, Matrix_DimSym);
1247        int nf = length(Omega), nz = INTEGER(Dim)[1];
1248    
1249        slot_dup(Matrix_DSym);
1250        slot_dup(Matrix_DIsqrtSym);
1251        slot_dup(Matrix_DimSym);
1252        slot_dup(Matrix_GpSym);
1253        slot_dup(Matrix_LiSym);
1254        slot_dup(Matrix_LpSym);
1255        slot_dup(Matrix_LxSym);
1256        slot_dup(Matrix_OmegaSym);
1257        slot_dup(Matrix_ParentSym);
1258        slot_dup(Matrix_bVarSym);
1259        slot_dup(Matrix_devianceSym);
1260        slot_dup(Matrix_devCompSym);
1261        slot_dup(Matrix_iSym);
1262        slot_dup(Matrix_ncSym);
1263        slot_dup(Matrix_statusSym);
1264        slot_dup(Matrix_pSym);
1265        slot_dup(Matrix_xSym);
1266        INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;
1267        SET_SLOT(ans, Matrix_XtXSym, allocMatrix(REALSXP, 1, 1));
1268        REAL(GET_SLOT(ans, Matrix_XtXSym))[0] = NA_REAL;
1269        SET_SLOT(ans, Matrix_RXXSym, allocMatrix(REALSXP, 1, 1));
1270        REAL(GET_SLOT(ans, Matrix_RXXSym))[0] = NA_REAL;
1271        SET_SLOT(ans, Matrix_ZtXSym, allocMatrix(REALSXP, nz, 1));
1272        SET_SLOT(ans, Matrix_RZXSym, allocMatrix(REALSXP, nz, 1));
1273        LOGICAL(GET_SLOT(ans, Matrix_statusSym))[0] = 0;
1274        UNPROTECT(1);
1275        return ans;
1276  }  }
1277    
1278    SEXP ssclme_to_lme(SEXP call, SEXP facs, SEXP x, SEXP model, SEXP REML,
1279                       SEXP rep, SEXP fitted, SEXP residuals)
1280    {
1281        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lme")));
1282    
1283        SET_SLOT(ans, install("call"), call);
1284        SET_SLOT(ans, install("facs"), facs);
1285        SET_SLOT(ans, Matrix_xSym, x);
1286        SET_SLOT(ans, install("model"), model);
1287        SET_SLOT(ans, install("REML"), REML);
1288        SET_SLOT(ans, install("rep"), rep);
1289        SET_SLOT(ans, install("fitted"), fitted);
1290        SET_SLOT(ans, install("residuals"), residuals);
1291        UNPROTECT(1);
1292        return ans;
1293    }

Legend:
Removed from v.100  
changed lines
  Added in v.164

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