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 21, Sat Mar 27 00:06:02 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 338  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 428  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 435  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 478  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 492  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 553  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 561  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 599  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 766  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 798  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 898  Line 922 
922      ssclme_invert(x);      ssclme_invert(x);
923      ryyinv = REAL(GET_SLOT(x, Matrix_RXXSym))[pp1*pp1 - 1];      ryyinv = REAL(GET_SLOT(x, Matrix_RXXSym))[pp1*pp1 - 1];
924      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
925          int nci = nc[i], Mi = (Gp[i+1] - Gp[i]), mi = Mi/nci;          int nci = nc[i], Mi = Gp[i+1] - Gp[i];
926          double *mm;          double *mm;
927    
928          SET_VECTOR_ELT(val, i, allocMatrix(REALSXP, nci, mi));          SET_VECTOR_ELT(val, i, allocMatrix(REALSXP, nci, Mi/nci));
929          mm = Memcpy(REAL(VECTOR_ELT(val, i)), b, Mi);          mm = Memcpy(REAL(VECTOR_ELT(val, i)), b, Mi);
930          b += nci * mi;          b += Mi;
931          for (j = 0; j < Mi; j++) mm[j] /= ryyinv;          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 1054  Line 1079 
1079                  *omgi = REAL(VECTOR_ELT(Omega, i)),                  *omgi = REAL(VECTOR_ELT(Omega, i)),
1080                  *tmp = Calloc(ncisq, double),                  *tmp = Calloc(ncisq, double),
1081                  diagj, one = 1.;                  diagj, one = 1.;
1082                /* FIXEME: Change this to use a factor and dsyrk */
1083                                  /* LD in omgi and L' in tmp */                                  /* LD in omgi and L' in tmp */
1084              memset(omgi, 0, sizeof(double) * ncisq);              memset(omgi, 0, sizeof(double) * ncisq);
1085              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
# Line 1165  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      }      }
# Line 1181  Line 1209 
1209      return R_NilValue;      return R_NilValue;
1210  }  }
1211    
1212  SEXP ssclme_asSscMatrix(SEXP x)  SEXP ssclme_gradient(SEXP x, SEXP REMLp)
1213  {  {
1214      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("sscMatrix")));      SEXP
1215      int *dims = INTEGER(GET_SLOT(val, Matrix_DimSym));          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      dims[0] = dims[1] = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];      p = pp1 - 1;
1237      SET_SLOT(val, Matrix_pSym, duplicate(GET_SLOT(x, Matrix_pSym)));      b = RZX + p * n;
1238      SET_SLOT(val, Matrix_iSym, duplicate(GET_SLOT(x, Matrix_iSym)));      ssclme_invert(x);
1239      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));      for (i = 0; i < nf; i++) {
1240      CHAR(STRING_ELT(GET_SLOT(x, Matrix_uploSym), 0))[0] = 'U';          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);      UNPROTECT(1);
1272      return val;      return ans;
1273  }  }
1274    
1275  SEXP ssclme_asTscMatrix(SEXP x)  SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats)
1276  {  {
1277      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("tscMatrix")));      SEXP val, b;
1278      int *dims = INTEGER(GET_SLOT(val, Matrix_DimSym));      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      dims[0] = dims[1] = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];      if (nf < 1)
1283      SET_SLOT(val, Matrix_pSym, duplicate(GET_SLOT(x, Matrix_LpSym)));          error("null factor list passed to ssclme_fitted");
1284      SET_SLOT(val, Matrix_iSym, duplicate(GET_SLOT(x, Matrix_LiSym)));      nobs = length(VECTOR_ELT(facs, 0));
1285      SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_LxSym)));      val = PROTECT(allocVector(REALSXP, nobs));
1286      CHAR(STRING_ELT(GET_SLOT(x, Matrix_uploSym), 0))[0] = 'U';      vv = REAL(val);
1287      CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0))[0] = 'U';      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);      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;      return val;
1351  }  }
1352    

Legend:
Removed from v.21  
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