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 10, Mon Mar 22 20:20:05 2004 UTC revision 100, Sat Apr 17 17:03:59 2004 UTC
# Line 52  Line 52 
52   *   *
53   * @return pointer to an integer R vector.   * @return pointer to an integer R vector.
54   */   */
55  static  
56  SEXP ctab_permute(SEXP ctab)  SEXP ctab_permute(SEXP ctab)
57  {  {
58      SEXP val, GpSl = GET_SLOT(ctab, Matrix_GpSym);      SEXP val, GpSl = GET_SLOT(ctab, Matrix_GpSym);
# Line 65  Line 65 
65          j,          j,
66          n = INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],          n = INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],
67          nf = length(GpSl) - 1,          nf = length(GpSl) - 1,
         nz = Ap[n],             /* number of non-zeros */  
68          pos;          pos;
69    
70      if (ctab_isNested(n, nf, 1, Ap, Ai, Gp))      if (ctab_isNested(n, nf, 1, Ap, Ai, Gp))
# Line 73  Line 72 
72      val =  allocVector(INTSXP, n);      val =  allocVector(INTSXP, n);
73      perm = INTEGER(val);      perm = INTEGER(val);
74      work = (int *) R_alloc(n, sizeof(int));      work = (int *) R_alloc(n, sizeof(int));
75      ssc_metis_order(n, nz, Ap, Ai, work, perm); /* perm gets inverse perm */      ssc_metis_order(n, Ap, Ai, work, perm);     /* perm gets inverse perm */
76      /* work now contains desired permutation but with groups scrambled */      /* work now contains desired permutation but with groups scrambled */
77    
78      /* copy work into perm preserving the order of the groups */      /* copy work into perm preserving the order of the groups */
# Line 189  Line 188 
188      }      }
189  }  }
190    
 static  
191  void ssclme_fill_LIp(int n, const int Parent[], int LIp[])  void ssclme_fill_LIp(int n, const int Parent[], int LIp[])
192  {  {
193      int *sz = Calloc(n, int), i;      int *sz = Calloc(n, int), i;
# Line 219  Line 217 
217  ssclme_create(SEXP facs, SEXP ncv, SEXP threshold)  ssclme_create(SEXP facs, SEXP ncv, SEXP threshold)
218  {  {
219      SEXP ctab, nms, ssc, tmp,      SEXP ctab, nms, ssc, tmp,
220          val = PROTECT(allocVector(VECSXP, 2));          val = PROTECT(allocVector(VECSXP, 2)),
221            dd = PROTECT(allocVector(INTSXP, 3));   /* dimensions of 3-D arrays */
222      int *Ai, *Ap, *Gp, *LIp, *Lp, *Parent,      int *Ai, *Ap, *Gp, *LIp, *Lp, *Parent,
223          *nc, Lnz, i, nf = length(facs), nzcol, pp1;          *nc, Lnz, i, nf = length(facs), nzcol, pp1,
224            *dims = INTEGER(dd);
225    
226      if (length(ncv) != (nf + 1))      if (length(ncv) != (nf + 1))
227          error("length of nc (%d) should be length of facs (%d) + 1",          error("length of nc (%d) should be length of facs (%d) + 1",
# Line 296  Line 296 
296      tmp = GET_SLOT(ssc, Matrix_bVarSym);      tmp = GET_SLOT(ssc, Matrix_bVarSym);
297      setAttrib(tmp, R_NamesSymbol, getAttrib(facs, R_NamesSymbol));      setAttrib(tmp, R_NamesSymbol, getAttrib(facs, R_NamesSymbol));
298      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
299          int ncol = Gp[i+1] - Gp[i], nrow = nc[i];          int nci = nc[i], mi = (Gp[i+1] - Gp[i])/nc[i];
300    
301          SET_VECTOR_ELT(tmp, i, allocMatrix(REALSXP, nrow, ncol));          dims[0] = dims[1] = nci;
302            dims[2] = mi;
303            SET_VECTOR_ELT(tmp, i, allocArray(REALSXP, dd));
304          memset(REAL(VECTOR_ELT(tmp, i)), 0,          memset(REAL(VECTOR_ELT(tmp, i)), 0,
305                 sizeof(double) * nrow * ncol);                 sizeof(double) * nci * nci * mi);
306      }      }
307      SET_SLOT(ssc, Matrix_LIpSym, allocVector(INTSXP, nzcol + 1));      SET_SLOT(ssc, Matrix_LIpSym, allocVector(INTSXP, nzcol + 1));
308      LIp = INTEGER(GET_SLOT(ssc, Matrix_LIpSym));      LIp = INTEGER(GET_SLOT(ssc, Matrix_LIpSym));
# Line 313  Line 315 
315          memset(REAL(GET_SLOT(ssc, Matrix_LIxSym)), 0,          memset(REAL(GET_SLOT(ssc, Matrix_LIxSym)), 0,
316                 sizeof(double) * Lnz);                 sizeof(double) * Lnz);
317      }      }
318      UNPROTECT(1);      UNPROTECT(2);
319      return val;      return val;
320  }  }
321    
# Line 334  Line 336 
336  }  }
337    
338  SEXP  SEXP
339  ssclme_update_mm(SEXP x, SEXP facs, SEXP mmats)  ssclme_transfer_dimnames(SEXP x, SEXP facs, SEXP mmats)
340  {  {
341        SEXP bVar = GET_SLOT(x, Matrix_bVarSym),
342            nms2 = PROTECT(allocVector(VECSXP, 2)),
343            nms3 = PROTECT(allocVector(VECSXP, 3));
344        int i, nf = length(mmats) - 1;
345        SEXP xcols = VECTOR_ELT(GetArrayDimnames(VECTOR_ELT(mmats, nf)), 1);
346    
347        for (i = 0; i < nf; i++) {
348            SEXP cnms = VECTOR_ELT(GetArrayDimnames(VECTOR_ELT(mmats, i)), 1);
349            SET_VECTOR_ELT(nms3, 0, cnms);
350            SET_VECTOR_ELT(nms3, 1, cnms);
351            SET_VECTOR_ELT(nms3, 2,
352                           getAttrib(VECTOR_ELT(facs, i), R_LevelsSymbol));
353            dimnamesgets(VECTOR_ELT(bVar, i), duplicate(nms3));
354        }
355        SET_VECTOR_ELT(nms2, 0, xcols);
356        SET_VECTOR_ELT(nms2, 1, xcols);
357        dimnamesgets(GET_SLOT(x, Matrix_XtXSym), nms2);
358        dimnamesgets(GET_SLOT(x, Matrix_RXXSym), nms2);
359        UNPROTECT(2);
360        return R_NilValue;
361    }
362    
363      SEXP      SEXP
364          bVar = GET_SLOT(x, Matrix_bVarSym);  ssclme_update_mm(SEXP x, SEXP facs, SEXP mmats)
365    {
366        SEXP bVar = GET_SLOT(x, Matrix_bVarSym);
367      int      int
368          *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),          *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
369          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
# Line 424  Line 450 
450                  if (ind < 0) error("logic error in ssclme_update_mm");                  if (ind < 0) error("logic error in ssclme_update_mm");
451                  if (Ncj || nck > 1) {                  if (Ncj || nck > 1) {
452                                  /* FIXME: run a loop to update */                                  /* FIXME: run a loop to update */
453                        error("code not yet written");
454                  } else {        /* update scalars directly */                  } else {        /* update scalars directly */
455                      Ax[ind] += Zj[fpji] * Zk[fpki];                      Ax[ind] += Zj[fpji] * Zk[fpki];
456                  }                  }
# Line 431  Line 458 
458          }          }
459      }      }
460      Free(Z);      Free(Z);
461        ssclme_transfer_dimnames(x, facs, mmats);
462      return R_NilValue;      return R_NilValue;
463  }  }
464    
465  SEXP ssclme_inflate_and_factor(SEXP lme)  SEXP ssclme_inflate_and_factor(SEXP x)
466  {  {
467      SEXP      SEXP
468          GpSlot = GET_SLOT(lme, Matrix_GpSym),          GpSlot = GET_SLOT(x, Matrix_GpSym),
469          Omega = GET_SLOT(lme, Matrix_OmegaSym);          Omega = GET_SLOT(x, Matrix_OmegaSym);
470      int n = INTEGER(GET_SLOT(lme, Matrix_DimSym))[1];      int n = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];
471      int      int
472          *Ai = INTEGER(GET_SLOT(lme, Matrix_iSym)),          *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
473          *Ap = INTEGER(GET_SLOT(lme, Matrix_pSym)),          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
474          *Flag = Calloc(n, int),          *Flag = Calloc(n, int),
475          *Gp = INTEGER(GpSlot),          *Gp = INTEGER(GpSlot),
476          *Lnz = Calloc(n, int),          *Lnz = Calloc(n, int),
477          *Pattern = Calloc(n, int),          *Pattern = Calloc(n, int),
478          *nc = INTEGER(GET_SLOT(lme, Matrix_ncSym)),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
479          j,          j,
480          nf = length(GpSlot) - 1;          nf = length(GpSlot) - 1;
481      double      double
482          *D = REAL(GET_SLOT(lme, Matrix_DSym)),          *D = REAL(GET_SLOT(x, Matrix_DSym)),
483          *DIsqrt = REAL(GET_SLOT(lme, Matrix_DIsqrtSym)),          *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
484          *Y = Calloc(n, double),          *Y = Calloc(n, double),
485          *xcp = Calloc(Ap[n], double);          *xcp = Calloc(Ap[n], double);
486    
487      Memcpy(xcp, REAL(GET_SLOT(lme, Matrix_xSym)), Ap[n]);      Memcpy(xcp, REAL(GET_SLOT(x, Matrix_xSym)), Ap[n]);
488      for (j = 0; j < nf; j++) {      for (j = 0; j < nf; j++) {
489          int  diag, i, ii, k, G2 = Gp[j + 1], ncj = nc[j];          int  diag, i, ii, k, G2 = Gp[j + 1], ncj = nc[j];
490          double *omgj = REAL(VECTOR_ELT(Omega, j));          double *omgj = REAL(VECTOR_ELT(Omega, j));
# Line 474  Line 502 
502          }          }
503      }      }
504      j = ldl_numeric(n, Ap, Ai, xcp,      j = ldl_numeric(n, Ap, Ai, xcp,
505                      INTEGER(GET_SLOT(lme, Matrix_LpSym)),                      INTEGER(GET_SLOT(x, Matrix_LpSym)),
506                      INTEGER(GET_SLOT(lme, Matrix_ParentSym)),                      INTEGER(GET_SLOT(x, Matrix_ParentSym)),
507                      Lnz, INTEGER(GET_SLOT(lme, Matrix_LiSym)),                      Lnz, INTEGER(GET_SLOT(x, Matrix_LiSym)),
508                      REAL(GET_SLOT(lme, Matrix_LxSym)),                      REAL(GET_SLOT(x, Matrix_LxSym)),
509                      D, Y, Pattern, Flag,                      D, Y, Pattern, Flag,
510                      (int *) NULL, (int *) NULL); /* P & Pinv */                      (int *) NULL, (int *) NULL); /* P & Pinv */
511      if (j != n)      if (j != n)
# Line 488  Line 516 
516      return R_NilValue;      return R_NilValue;
517  }  }
518    
519  SEXP ssclme_factor(SEXP lme)  SEXP ssclme_factor(SEXP x)
520  {  {
521      int *status = LOGICAL(GET_SLOT(lme, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
522    
523      if (!status[0]) {      if (!status[0]) {
524          SEXP          SEXP
525              GpSlot = GET_SLOT(lme, Matrix_GpSym),              GpSlot = GET_SLOT(x, Matrix_GpSym),
526              Omega = GET_SLOT(lme, Matrix_OmegaSym);              Omega = GET_SLOT(x, Matrix_OmegaSym);
527          int          int
528              *Gp = INTEGER(GpSlot),              *Gp = INTEGER(GpSlot),
529              *Li = INTEGER(GET_SLOT(lme, Matrix_LiSym)),              *Li = INTEGER(GET_SLOT(x, Matrix_LiSym)),
530              *Lp = INTEGER(GET_SLOT(lme, Matrix_LpSym)),              *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
531              *nc = INTEGER(GET_SLOT(lme, Matrix_ncSym)),              *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
532              i,              i,
533              n = INTEGER(GET_SLOT(lme, Matrix_DimSym))[1],              n = INTEGER(GET_SLOT(x, Matrix_DimSym))[1],
534              nf = length(GpSlot) - 1,              nf = length(GpSlot) - 1,
535              nobs = nc[nf + 1],              nobs = nc[nf + 1],
536              nreml = nobs + 1 - nc[nf],              nreml = nobs + 1 - nc[nf],
537              pp1 = nc[nf],              pp1 = nc[nf],
538              pp2 = pp1 + 1;              pp2 = pp1 + 1;
539          double          double
540              *D = REAL(GET_SLOT(lme, Matrix_DSym)),              *D = REAL(GET_SLOT(x, Matrix_DSym)),
541              *DIsqrt = REAL(GET_SLOT(lme, Matrix_DIsqrtSym)),              *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
542              *Lx = REAL(GET_SLOT(lme, Matrix_LxSym)),              *Lx = REAL(GET_SLOT(x, Matrix_LxSym)),
543              *RXX = REAL(GET_SLOT(lme, Matrix_RXXSym)),              *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
544              *RZX = REAL(GET_SLOT(lme, Matrix_RZXSym)),              *RZX = REAL(GET_SLOT(x, Matrix_RZXSym)),
545              *dcmp = REAL(getAttrib(lme, Matrix_devCompSym)),              *dcmp = REAL(getAttrib(x, Matrix_devCompSym)),
546              *deviance = REAL(getAttrib(lme, Matrix_devianceSym)),              *deviance = REAL(getAttrib(x, Matrix_devianceSym)),
547              minus1 = -1.,              minus1 = -1.,
548              one = 1.;              one = 1.;
549    
550          ssclme_inflate_and_factor(lme);          ssclme_inflate_and_factor(x);
551                                  /* Accumulate logdet of ZtZ+W */                                  /* Accumulate logdet of ZtZ+W */
552          dcmp[0] = dcmp[1] = dcmp[2] = dcmp[3] = 0.;          dcmp[0] = dcmp[1] = dcmp[2] = dcmp[3] = 0.;
553          for (i = 0; i < n; i++) dcmp[0] += log(D[i]);          for (i = 0; i < n; i++) dcmp[0] += log(D[i]);
# Line 549  Line 577 
577              }              }
578          }          }
579                                  /* ldl_lsolve on Z'X */                                  /* ldl_lsolve on Z'X */
580          Memcpy(RZX, REAL(GET_SLOT(lme, Matrix_ZtXSym)), n * pp1);          Memcpy(RZX, REAL(GET_SLOT(x, Matrix_ZtXSym)), n * pp1);
581          for (i = 0; i < pp1; i++) {          for (i = 0; i < pp1; i++) {
582              int j;              int j;
583              double *RZXi = RZX + i * n;              double *RZXi = RZX + i * n;
# Line 557  Line 585 
585              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
586          }          }
587                                  /* downdate and factor X'X */                                  /* downdate and factor X'X */
588          Memcpy(RXX, REAL(GET_SLOT(lme, Matrix_XtXSym)), pp1 * pp1);          Memcpy(RXX, REAL(GET_SLOT(x, Matrix_XtXSym)), pp1 * pp1);
589          F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,          F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,
590                          RZX, &n, &one, RXX, &pp1);                          RZX, &n, &one, RXX, &pp1);
591          F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);          F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);
# Line 595  Line 623 
623   * @return R_NilValue (x is updated in place)   * @return R_NilValue (x is updated in place)
624    
625   */   */
626    static
627  SEXP ldl_inverse(SEXP x)  SEXP ldl_inverse(SEXP x)
628  {  {
629      SEXP      SEXP
# Line 762  Line 791 
791      }      }
792      return R_NilValue;      return R_NilValue;
793  }  }
794    SEXP ssclme_invert(SEXP x)
 SEXP ssclme_invert(SEXP lme)  
795  {  {
796      int *status = LOGICAL(GET_SLOT(lme, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
797      if (!status[0]) ssclme_factor(lme);      if (!status[0]) ssclme_factor(x);
798      if (!status[1]) {      if (!status[1]) {
799          SEXP          SEXP
800              RZXsl = GET_SLOT(lme, Matrix_RZXSym);              RZXsl = GET_SLOT(x, Matrix_RZXSym);
801          int          int
802              *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),              *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
803              *Li = INTEGER(GET_SLOT(lme, Matrix_LiSym)),              *Li = INTEGER(GET_SLOT(x, Matrix_LiSym)),
804              *Lp = INTEGER(GET_SLOT(lme, Matrix_LpSym)),              *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
805              i,              i,
806              n = dims[0],              n = dims[0],
807              pp1 = dims[1];              pp1 = dims[1];
808          double          double
809              *DIsqrt = REAL(GET_SLOT(lme, Matrix_DIsqrtSym)),              *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
810              *Lx = REAL(GET_SLOT(lme, Matrix_LxSym)),              *Lx = REAL(GET_SLOT(x, Matrix_LxSym)),
811              *RXX = REAL(GET_SLOT(lme, Matrix_RXXSym)),              *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
812              *RZX = REAL(RZXsl),              *RZX = REAL(RZXsl),
813              one = 1.;              one = 1.;
814    
# Line 794  Line 822 
822              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
823              ldl_ltsolve(n, RZXi, Lp, Li, Lx);              ldl_ltsolve(n, RZXi, Lp, Li, Lx);
824          }          }
825          ldl_inverse(lme);          ldl_inverse(x);
826          status[1] = 1;          status[1] = 1;
827      }      }
828      return R_NilValue;      return R_NilValue;
# Line 876  Line 904 
904   */   */
905  SEXP ssclme_ranef(SEXP x)  SEXP ssclme_ranef(SEXP x)
906  {  {
907      SEXP RZXsl = GET_SLOT(x, Matrix_RZXSym);      SEXP RZXsl = GET_SLOT(x, Matrix_RZXSym),
908            GpSl = GET_SLOT(x, Matrix_GpSym);
909      int *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),      int *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
910          j,          *Gp = INTEGER(GpSl),
911            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
912            i, j,
913          n = dims[0],          n = dims[0],
914            nf = length(GpSl) - 1,
915          pp1 = dims[1],          pp1 = dims[1],
916          p = pp1 - 1;          p = pp1 - 1;
917      SEXP val = PROTECT(allocVector(REALSXP, n));      SEXP val = PROTECT(allocVector(VECSXP, nf));
918      double      double
919          *RZX = REAL(RZXsl),          *b = REAL(RZXsl) + n * p,
         *b = REAL(val),  
920          ryyinv;         /* ryy-inverse */          ryyinv;         /* ryy-inverse */
921    
922      ssclme_invert(x);      ssclme_invert(x);
     Memcpy(b, RZX + p * n, n);  
923      ryyinv = REAL(GET_SLOT(x, Matrix_RXXSym))[pp1*pp1 - 1];      ryyinv = REAL(GET_SLOT(x, Matrix_RXXSym))[pp1*pp1 - 1];
924      for (j = 0; j < n; j++) b[j] /= ryyinv;      for (i = 0; i < nf; i++) {
925            int nci = nc[i], Mi = Gp[i+1] - Gp[i];
926            double *mm;
927    
928            SET_VECTOR_ELT(val, i, allocMatrix(REALSXP, nci, Mi/nci));
929            mm = Memcpy(REAL(VECTOR_ELT(val, i)), b, Mi);
930            b += Mi;
931            for (j = 0; j < Mi; j++) mm[j] /= ryyinv;
932        }
933      UNPROTECT(1);      UNPROTECT(1);
934      return val;      return val;
935  }  }
936    
937  /**  /**
938   * Extract the ML or REML conditional estimate of sigma   * Extract the ML or REML conditional estimate of sigma
939   *   *
# Line 957  Line 996 
996  }  }
997    
998  /**  /**
999     * Extract the upper triangles of the Omega matrices in the unconstrained
1000     * parameterization.
1001     * (These are not in any sense "coefficients" but the extractor is
1002     * called coef for historical reasons.)
1003     *
1004     * @param x pointer to an ssclme object
1005     *
1006     * @return numeric vector of the values in the upper triangles of the
1007     * Omega matrices
1008     */
1009    SEXP ssclme_coefUnc(SEXP x)
1010    {
1011        SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1012        int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1013            i, nf = length(Omega), vind;
1014        SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
1015        double *vv = REAL(val);
1016    
1017        vind = 0;
1018        for (i = 0; i < nf; i++) {
1019            int  nci = nc[i];
1020            if (nci == 1) {
1021                vv[vind++] = log(REAL(VECTOR_ELT(Omega, i))[0]);
1022            } else {
1023                int j, k, ncip1 = nci + 1, ncisq = nci * nci;
1024                double *tmp = Memcpy(Calloc(ncisq, double),
1025                                     REAL(VECTOR_ELT(Omega, i)), ncisq);
1026                F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
1027                if (j)              /* should never happen */
1028                    error("DPOTRF returned error code %d", j);
1029                for (j = 0; j < nci; j++) {
1030                    double diagj = tmp[j * ncip1];
1031                    vv[vind++] = 2. * log(diagj);
1032                    for (k = j + 1; k < nci; k++) {
1033                        tmp[j + k * nci] /= diagj;
1034                    }
1035                }
1036                for (j = 0; j < nci; j++) {
1037                    for (k = j + 1; k < nci; k++) {
1038                        vv[vind++] = tmp[j + k * nci];
1039                    }
1040                }
1041                Free(tmp);
1042            }
1043        }
1044        UNPROTECT(1);
1045        return val;
1046    }
1047    
1048    /**
1049     * Assign the upper triangles of the Omega matrices in the unconstrained
1050     * parameterization.
1051     *
1052     * @param x pointer to an ssclme object
1053     * @param coef pointer to an numeric vector of appropriate length
1054     *
1055     * @return R_NilValue
1056     */
1057    SEXP ssclme_coefGetsUnc(SEXP x, SEXP coef)
1058    {
1059        SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1060        int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1061            cind, i, nf = length(Omega),
1062            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
1063        double *cc = REAL(coef);
1064    
1065        if (length(coef) != coef_length(nf, nc) || !isReal(coef))
1066            error("coef must be a numeric vector of length %d",
1067                  coef_length(nf, nc));
1068        cind = 0;
1069        for (i = 0; i < nf; i++) {
1070            int nci = nc[i];
1071            if (nci == 1) {
1072                REAL(VECTOR_ELT(Omega, i))[0] = exp(cc[cind++]);
1073            } else {
1074                int odind = cind + nci, /* off-diagonal index */
1075                    j, k,
1076                    ncip1 = nci + 1,
1077                    ncisq = nci * nci;
1078                double
1079                    *omgi = REAL(VECTOR_ELT(Omega, i)),
1080                    *tmp = Calloc(ncisq, double),
1081                    diagj, one = 1.;
1082                /* FIXEME: Change this to use a factor and dsyrk */
1083                                    /* LD in omgi and L' in tmp */
1084                memset(omgi, 0, sizeof(double) * ncisq);
1085                for (j = 0; j < nci; j++) {
1086                    omgi[j * ncip1] = diagj = exp(cc[cind++]);
1087                    for (k = j + 1; k < nci; k++) {
1088                        omgi[j*nci + k] = diagj * (tmp[k*nci + j] = cc[odind++]);
1089                    }
1090                }
1091                F77_CALL(dtrmm)("R", "U", "N", "U", &nci, &nci, &one,
1092                                tmp, &nci, omgi, &nci);
1093                Free(tmp);
1094                cind = odind;
1095            }
1096        }
1097        status[0] = status[1] = 0;
1098        return x;
1099    }
1100    
1101    /**
1102   * Assign the upper triangles of the Omega matrices.   * Assign the upper triangles of the Omega matrices.
1103   * (These are not in any sense "coefficients" but are   * (These are not in any sense "coefficients" but are
1104   * called coef for historical reasons.)   * called coef for historical reasons.)
# Line 988  Line 1130 
1130          }          }
1131      }      }
1132      status[0] = status[1] = 0;      status[0] = status[1] = 0;
1133      return R_NilValue;      return x;
1134  }  }
1135    
1136  SEXP ssclme_EMstepsGets(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)  SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)
1137  {  {
1138      SEXP      SEXP
1139          Omega = GET_SLOT(x, Matrix_OmegaSym),          Omega = GET_SLOT(x, Matrix_OmegaSym),
# Line 1049  Line 1191 
1191                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,
1192                              &one, vali, &nci);                              &one, vali, &nci);
1193              if (REML) {              if (REML) {
1194                  int mp = mi * p;                  int j;
1195                  F77_CALL(dsyrk)("U", "N", &nci, &mp,                  for (j = 0; j < p; j++) {
1196                                  &alpha, RZX + Gp[i], &nci,                      F77_CALL(dsyrk)("U", "N", &nci, &mi,
1197                                    &alpha, RZX + Gp[i] + j*n, &nci,
1198                                  &one, vali, &nci);                                  &one, vali, &nci);
1199              }              }
1200                }
1201              F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);              F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);
1202              if (info) error("DPOTRF returned error code %d", info);              if (info) error("DPOTRF returned error code %d", info);
1203              F77_CALL(dpotri)("U", &nci, vali, &nci, &info);              F77_CALL(dpotri)("U", &nci, vali, &nci, &info);
1204              if (info) error("DPOTRF returned error code %d", info);              if (info) error("DPOTRI returned error code %d", info);
1205          }          }
1206          status[0] = status[1] = 0;          status[0] = status[1] = 0;
1207      }      }
1208      ssclme_factor(x);      ssclme_factor(x);
1209      return R_NilValue;      return R_NilValue;
1210  }  }
1211    
1212    SEXP ssclme_gradient(SEXP x, SEXP REMLp)
1213    {
1214        SEXP
1215            RZXsl = GET_SLOT(x, Matrix_RZXSym),
1216            ans = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym))),
1217            ncsl = GET_SLOT(x, Matrix_ncSym),
1218            bVar = GET_SLOT(x, Matrix_bVarSym);
1219        int
1220            *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1221            *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1222            *nc = INTEGER(ncsl),
1223            REML = asLogical(REMLp),
1224            i, info,
1225            n = dims[0],
1226            nf = length(ncsl) - 2,
1227            nobs = nc[nf + 1],
1228            p,
1229            pp1 = dims[1];
1230        double
1231            *RZX = REAL(RZXsl),
1232            *b,
1233            alpha,
1234            one = 1.;
1235    
1236        p = pp1 - 1;
1237        b = RZX + p * n;
1238        ssclme_invert(x);
1239        for (i = 0; i < nf; i++) {
1240            int ki = Gp[i+1] - Gp[i],
1241                nci = nc[i],
1242                mi = ki/nci;
1243            double
1244                *vali = REAL(VECTOR_ELT(ans, i));
1245    
1246            F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);
1247            if (info)
1248                error("DPOTRF returned error code %d for component %d of Omega",
1249                      info, i + 1);
1250            F77_CALL(dpotri)("U", &nci, vali, &nci, &info);
1251            if (info)
1252                error("DPOTRI returned error code %d for component %d of Omega",
1253                      info, i + 1);
1254            alpha = (double) -mi;
1255            F77_CALL(dsyrk)("U", "N", &nci, &ki,
1256                            &one, REAL(VECTOR_ELT(bVar, i)), &nci,
1257                            &alpha, vali, &nci);
1258            alpha = ((double)(REML?(nobs-p):nobs));
1259            F77_CALL(dsyrk)("U", "N", &nci, &mi,
1260                            &alpha, b + Gp[i], &nci,
1261                            &one, vali, &nci);
1262            if (REML) {
1263                int j;
1264                for (j = 0; j < p; j++) {
1265                    F77_CALL(dsyrk)("U", "N", &nci, &mi,
1266                                    &one, RZX + Gp[i] + j*n, &nci,
1267                                    &one, vali, &nci);
1268                }
1269            }
1270        }
1271        UNPROTECT(1);
1272        return ans;
1273    }
1274    
1275    SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats)
1276    {
1277        SEXP val, b;
1278        int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1279            i, ione = 1, nf = length(facs), nobs, p;
1280        double *vv, one = 1.0, zero = 0.0;
1281    
1282        if (nf < 1)
1283            error("null factor list passed to ssclme_fitted");
1284        nobs = length(VECTOR_ELT(facs, 0));
1285        val = PROTECT(allocVector(REALSXP, nobs));
1286        vv = REAL(val);
1287        p = nc[nf] - 1;
1288        if (p > 0) {
1289            F77_CALL(dgemm)("N", "N", &nobs, &ione, &p, &one,
1290                            REAL(VECTOR_ELT(mmats, nf)), &nobs,
1291                            REAL(PROTECT(ssclme_fixef(x))), &p,
1292                            &zero, vv, &nobs);
1293            UNPROTECT(1);
1294        } else {
1295            memset(vv, 0, sizeof(double) * nobs);
1296        }
1297        b = PROTECT(ssclme_ranef(x));
1298        for (i = 0; i < nf; i++) {
1299            int *ff = INTEGER(VECTOR_ELT(facs, i)), j, nci = nc[i];
1300            double *bb = REAL(VECTOR_ELT(b, i)),
1301                *mm = REAL(VECTOR_ELT(mmats, i));
1302            for (j = 0; j < nobs; ) {
1303                int nn = 1, lev = ff[j];
1304                /* check for adjacent rows with same factor level */
1305                while (ff[j + nn] == lev) nn++;
1306                F77_CALL(dgemm)("N", "N", &nn, &ione, &nci,
1307                                &one, mm + j, &nobs,
1308                                bb + (lev - 1) * nci, &nci,
1309                                &one, vv + j, &nobs);
1310                j += nn;
1311            }
1312        }
1313        UNPROTECT(2);
1314        return val;
1315    }
1316    
1317    SEXP ssclme_variances(SEXP x, SEXP REML)
1318    {
1319        SEXP Omg = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym))),
1320            val = PROTECT(allocVector(VECSXP, 2));
1321        int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1322            i, nf = length(Omg);
1323        double sigmasq;
1324    
1325    
1326        SET_VECTOR_ELT(val, 0, Omg);
1327        SET_VECTOR_ELT(val, 1, ssclme_sigma(x, REML));
1328        sigmasq = REAL(VECTOR_ELT(val, 1))[0];
1329        sigmasq = (sigmasq) * (sigmasq);
1330        for (i = 0; i < nf; i++) {
1331            double *mm = REAL(VECTOR_ELT(Omg, i));
1332            int j, k, nci = nc[i], ncip1 = nci+1;
1333    
1334            F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);
1335            if (j)                  /* shouldn't happen */
1336                error("DPOTRF returned error code %d on Omega[%d]",
1337                      j, i + 1);
1338            F77_CALL(dpotri)("U", &nci, mm, &nci, &j);
1339            if (j)                  /* shouldn't happen */
1340                error("DTRTRI returned error code %d on Omega[%d]",
1341                      j, i + 1);
1342            for (j = 0; j < nci; j++) {
1343                mm[j * ncip1] *= sigmasq;
1344                for (k = 0; k < j; k++) {
1345                    mm[j + k * nci] = (mm[k + j * nci] *= sigmasq);
1346                }
1347            }
1348        }
1349        UNPROTECT(2);
1350        return val;
1351    }
1352    

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

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