# SCM Repository

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

# Diff of /pkg/src/ssclme.c

revision 103, Sun Apr 18 18:31:43 2004 UTC revision 104, Sun Apr 18 19:21:57 2004 UTC
# Line 868  Line 868
868
869  /**  /**
870   * Extract the conditional estimates of the fixed effects   * Extract the conditional estimates of the fixed effects
871   *   *
872   * @param x Pointer to an ssclme object   * @param x Pointer to an ssclme object
873   *   *
# Line 896  Line 895
895
896  /**  /**
897   * 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.
898   *   *
899   * @param x Pointer to an ssclme object   * @param x Pointer to an ssclme object
900   *   *
# Line 966  Line 963
963
964  /**  /**
965   * Extract the upper triangles of the Omega matrices.   * Extract the upper triangles of the Omega matrices.
966   * (These are not in any sense "coefficients" but the extractor is   * (These aren't "coefficients" but the extractor is
967   * called coef for historical reasons.)   * called coef for historical reasons.)
968   *   *
969   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
# Line 984  Line 981
981
982      vind = 0;      vind = 0;
983      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
984          int j, k, nci = nc[i];          int nci = nc[i];
985            if (nci == 1) {
986                vv[vind++] = REAL(VECTOR_ELT(Omega, i))[0];
987            } else {
988                int j, k, odind = vind + nci, ncip1 = nci + 1;
989          double *omgi = REAL(VECTOR_ELT(Omega, i));          double *omgi = REAL(VECTOR_ELT(Omega, i));
990
991          for (j = 0; j < nci; j++) {          for (j = 0; j < nci; j++) {
992              for (k = 0; k <= j; k++) {                  vv[vind++] = omgi[j * ncip1];
993                  vv[vind++] = omgi[j*nci + k];                  for (k = j + 1; k < nci; k++) {
994                        vv[odind++] = omgi[k*nci + j];
995              }              }
996          }          }
997                vind = odind;
998            }
999      }      }
1000      UNPROTECT(1);      UNPROTECT(1);
1001      return val;      return val;
1002  }  }
1003
1004  /**  /**
1005   * Extract the upper triangles of the Omega matrices in the unconstrained   * Extract the unconstrained parameters that determine the
1006   * parameterization.   * Omega matrices. (Called coef for historical reasons.)
* (These are not in any sense "coefficients" but the extractor is
* called coef for historical reasons.)
1007   *   *
1008   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1009   *   *
1010   * @return numeric vector of the values in the upper triangles of the   * @return numeric vector of unconstrained parameters that determine the
1011   * Omega matrices   * Omega matrices
1012   */   */
1013  SEXP ssclme_coefUnc(SEXP x)  SEXP ssclme_coefUnc(SEXP x)
# Line 1026  Line 1029
1029                                   REAL(VECTOR_ELT(Omega, i)), ncisq);                                   REAL(VECTOR_ELT(Omega, i)), ncisq);
1030              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
1031              if (j)              /* should never happen */              if (j)              /* should never happen */
1032                  error("DPOTRF returned error code %d", j);                  error("DPOTRF returned error code %d on Omega[[%d]]",
1033                          j, i+1);
1034              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
1035                  double diagj = tmp[j * ncip1];                  double diagj = tmp[j * ncip1];
1036                  vv[vind++] = 2. * log(diagj);                  vv[vind++] = 2. * log(diagj);
# Line 1047  Line 1051
1051  }  }
1052
1053  /**  /**
1054   * Assign the upper triangles of the Omega matrices in the unconstrained   * Assign the Omega matrices from the unconstrained parameterization.
* parameterization.
1055   *   *
1056   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1057   * @param coef pointer to an numeric vector of appropriate length   * @param coef pointer to an numeric vector of appropriate length
# Line 1079  Line 1082
1082              double              double
1083                  *omgi = REAL(VECTOR_ELT(Omega, i)),                  *omgi = REAL(VECTOR_ELT(Omega, i)),
1084                  *tmp = Calloc(ncisq, double),                  *tmp = Calloc(ncisq, double),
1085                  diagj, one = 1.;                  diagj, one = 1., zero = 0.;
1086              /* FIXEME: Change this to use a factor and dsyrk */
/* LD in omgi and L' in tmp */
1087              memset(omgi, 0, sizeof(double) * ncisq);              memset(omgi, 0, sizeof(double) * ncisq);
1088              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
1089                  omgi[j * ncip1] = diagj = exp(cc[cind++]);                  tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);
1090                  for (k = j + 1; k < nci; k++) {                  for (k = j + 1; k < nci; k++) {
1091                      omgi[j*nci + k] = diagj * (tmp[k*nci + j] = cc[odind++]);                      tmp[k*nci + j] = cc[odind++] * diagj;
1092                  }                  }
1093              }              }
1094              F77_CALL(dtrmm)("R", "U", "N", "U", &nci, &nci, &one,              F77_CALL(dsyrk)("U", "T", &nci, &nci, &one,
1095                              tmp, &nci, omgi, &nci);                              tmp, &nci, &zero, omgi, &nci);
1096              Free(tmp);              Free(tmp);
1097              cind = odind;              cind = odind;
1098          }          }
# Line 1101  Line 1103
1103
1104  /**  /**
1105   * Assign the upper triangles of the Omega matrices.   * Assign the upper triangles of the Omega matrices.
1106   * (These are not in any sense "coefficients" but are   * (Called coef for historical reasons.)
* called coef for historical reasons.)
1107   *   *
1108   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1109   * @param coef pointer to an numeric vector of appropriate length   * @param coef pointer to an numeric vector of appropriate length
# Line 1122  Line 1123
1123                coef_length(nf, nc));                coef_length(nf, nc));
1124      cind = 0;      cind = 0;
1125      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1126          int j, k, nci = nc[i];          int nci = nc[i];
1127            if (nci == 1) {
1128                REAL(VECTOR_ELT(Omega, i))[0] = cc[cind++];
1129            } else {
1130                int j, k, odind = cind + nci, ncip1 = nci + 1;
1131          double *omgi = REAL(VECTOR_ELT(Omega, i));          double *omgi = REAL(VECTOR_ELT(Omega, i));
1132
1133          for (j = 0; j < nci; j++) {          for (j = 0; j < nci; j++) {
1134              for (k = 0; k <= j; k++) {                  omgi[j * ncip1] = cc[cind++];
1135                  omgi[j*nci + k] = cc[cind++];                  for (k = j + 1; k < nci; k++) {
1136                        omgi[k*nci + j] = cc[odind++];
1137              }              }
1138          }          }
1139                cind = odind;
1140            }
1141      }      }
1142      status[0] = status[1] = 0;      status[0] = status[1] = 0;
1143      return x;      return x;
# Line 1215  Line 1224
1224      SEXP      SEXP
1225          Omega = GET_SLOT(x, Matrix_OmegaSym),          Omega = GET_SLOT(x, Matrix_OmegaSym),
1226          RZXsl = GET_SLOT(x, Matrix_RZXSym),          RZXsl = GET_SLOT(x, Matrix_RZXSym),
1227          ans = PROTECT(duplicate(Omega)),          ans,
1228          ncsl = GET_SLOT(x, Matrix_ncSym),          ncsl = GET_SLOT(x, Matrix_ncSym),
1229          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
1230      int      int
# Line 1223  Line 1232
1232          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1233          *nc = INTEGER(ncsl),          *nc = INTEGER(ncsl),
1234          REML = asLogical(REMLp),          REML = asLogical(REMLp),
1235          i, info,          cind, i, n = dims[0],
n = dims[0],
1236          nf = length(ncsl) - 2,          nf = length(ncsl) - 2,
1237          nobs = nc[nf + 1],          nobs, odind, p, pp1 = dims[1],
p,
pp1 = dims[1],
1238          uncst = asLogical(Uncp);          uncst = asLogical(Uncp);
1239      double      double
1240          *RZX = REAL(RZXsl),          *RZX = REAL(RZXsl),
# Line 1236  Line 1242
1242          alpha,          alpha,
1243          one = 1.;          one = 1.;
1244
1245        nobs = nc[nf + 1];
1246      p = pp1 - 1;      p = pp1 - 1;
1247      b = RZX + p * n;      b = RZX + p * n;
1248        ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
1249      ssclme_invert(x);      ssclme_invert(x);
1250        cind = 0;
1251      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1252          int ki = Gp[i+1] - Gp[i],          int j, ki = Gp[i+1] - Gp[i],
1253              nci = nc[i],              nci = nc[i],
1254              mi = ki/nci;              mi = ki/nci;
1255          double          double
1256              *vali = REAL(VECTOR_ELT(ans, i));              *tmp = Memcpy(Calloc(nci * nci, double),
1257                              REAL(VECTOR_ELT(Omega, i)), nci * nci);
1258
1259          F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);          F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
1260          if (info)          if (j)
1261              error("DPOTRF returned error code %d for component %d of Omega",              error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1262                    info, i + 1);          F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);
1263          F77_CALL(dpotri)("U", &nci, vali, &nci, &info);          if (j)
1264          if (info)              error("DPOTRI gave error code %d on Omega[[%d]]", j, i + 1);
error("DPOTRI returned error code %d for component %d of Omega",
info, i + 1);
1265          alpha = (double) -mi;          alpha = (double) -mi;
1266          F77_CALL(dsyrk)("U", "N", &nci, &ki,          F77_CALL(dsyrk)("U", "N", &nci, &ki,
1267                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,
1268                          &alpha, vali, &nci);                          &alpha, tmp, &nci);
1269          alpha = ((double)(REML?(nobs-p):nobs));          alpha = ((double)(REML?(nobs-p):nobs));
1270          F77_CALL(dsyrk)("U", "N", &nci, &mi,          F77_CALL(dsyrk)("U", "N", &nci, &mi,
1271                          &alpha, b + Gp[i], &nci,                          &alpha, b + Gp[i], &nci,
1272                          &one, vali, &nci);                          &one, tmp, &nci);
1273          if (REML) {          if (REML) {
int j;
1274              for (j = 0; j < p; j++) {              for (j = 0; j < p; j++) {
1275                  F77_CALL(dsyrk)("U", "N", &nci, &mi,                  F77_CALL(dsyrk)("U", "N", &nci, &mi,
1276                                  &one, RZX + Gp[i] + j*n, &nci,                                  &one, RZX + Gp[i] + j*n, &nci,
1277                                  &one, vali, &nci);                                  &one, tmp, &nci);
1278              }              }
1279          }          }
if (uncst) {
1280              if (nci == 1) {              if (nci == 1) {
1281                  *vali *= *REAL(VECTOR_ELT(Omega, i));              REAL(ans)[cind++] = *tmp *
1282                    (uncst ? *REAL(VECTOR_ELT(Omega, i)) : 1.);
1283              } else {              } else {
1284                if (uncst) {
1285                  error("Code not written yet");                  error("Code not written yet");
1286                } else {
1287                    int k, ncip1 = nci + 1, odind = cind + nci;
1288                    for (j = 0; j < nci; j++) {
1289                        REAL(ans)[cind++] = tmp[j * ncip1];
1290                        for (k = j + 1; k < nci; k++) {
1291                            REAL(ans)[odind++] = tmp[k*nci + j] * 2.;
1292              }              }
1293          }          }
1294                    cind = odind;
1295                }
1296            }
1297            Free(tmp);
1298      }      }
1299      UNPROTECT(1);      UNPROTECT(1);
1300      return ans;      return ans;

Legend:
 Removed from v.103 changed lines Added in v.104