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 106, Mon Apr 19 00:48:02 2004 UTC
# Line 252  Line 252 
252      SET_SLOT(ssc, Matrix_RXXSym, allocMatrix(REALSXP, pp1, pp1));      SET_SLOT(ssc, Matrix_RXXSym, allocMatrix(REALSXP, pp1, pp1));
253      SET_SLOT(ssc, Matrix_ZtXSym, allocMatrix(REALSXP, nzcol, pp1));      SET_SLOT(ssc, Matrix_ZtXSym, allocMatrix(REALSXP, nzcol, pp1));
254      SET_SLOT(ssc, Matrix_RZXSym, allocMatrix(REALSXP, nzcol, pp1));      SET_SLOT(ssc, Matrix_RZXSym, allocMatrix(REALSXP, nzcol, pp1));
255          /* Zero the symmetric matrices (for cosmetic reasons only). */                                  /* Zero symmetric matrices (cosmetic) */
256      memset(REAL(GET_SLOT(ssc, Matrix_XtXSym)), 0,      memset(REAL(GET_SLOT(ssc, Matrix_XtXSym)), 0,
257             sizeof(double) * pp1 * pp1);             sizeof(double) * pp1 * pp1);
258      memset(REAL(GET_SLOT(ssc, Matrix_RXXSym)), 0,      memset(REAL(GET_SLOT(ssc, Matrix_RXXSym)), 0,
# Line 791  Line 791 
791      }      }
792      return R_NilValue;      return R_NilValue;
793  }  }
794    
795  SEXP ssclme_invert(SEXP x)  SEXP ssclme_invert(SEXP x)
796  {  {
797      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
# Line 867  Line 868 
868    
869  /**  /**
870   * Extract the conditional estimates of the fixed effects   * Extract the conditional estimates of the fixed effects
  * FIXME: Add names  
871   *   *
872   * @param x Pointer to an ssclme object   * @param x Pointer to an ssclme object
873   *   *
# Line 895  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 965  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 983  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);
# Line 996  Line 1002 
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 1025  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 1046  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 1078  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 1100  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 1121  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 1209  Line 1219 
1219      return R_NilValue;      return R_NilValue;
1220  }  }
1221    
1222  SEXP ssclme_gradient(SEXP x, SEXP REMLp)  SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)
1223  {  {
1224      SEXP      SEXP
1225            Omega = GET_SLOT(x, Matrix_OmegaSym),
1226          RZXsl = GET_SLOT(x, Matrix_RZXSym),          RZXsl = GET_SLOT(x, Matrix_RZXSym),
         ans = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym))),  
1227          ncsl = GET_SLOT(x, Matrix_ncSym),          ncsl = GET_SLOT(x, Matrix_ncSym),
1228          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
1229      int      int
# Line 1221  Line 1231 
1231          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1232          *nc = INTEGER(ncsl),          *nc = INTEGER(ncsl),
1233          REML = asLogical(REMLp),          REML = asLogical(REMLp),
1234          i, info,          cind, i, n = dims[0],
1235          n = dims[0],          nf = length(Omega),
1236          nf = length(ncsl) - 2,          nobs, odind, p, pp1 = dims[1],
1237          nobs = nc[nf + 1],          uncst = asLogical(Uncp);
         p,  
         pp1 = dims[1];  
1238      double      double
1239          *RZX = REAL(RZXsl),          *RZX = REAL(RZXsl),
1240          *b,          *b,
1241          alpha,          alpha,
1242          one = 1.;          one = 1.;
1243        SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
1244    
1245        nobs = nc[nf + 1];
1246      p = pp1 - 1;      p = pp1 - 1;
1247      b = RZX + p * n;      b = RZX + p * n;
1248      ssclme_invert(x);      ssclme_invert(x);
1249        cind = 0;
1250      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1251          int ki = Gp[i+1] - Gp[i],          int j, ki = Gp[i+1] - Gp[i],
1252              nci = nc[i],              nci = nc[i], ncip1 = nci + 1, ncisq = nci * nci,
1253              mi = ki/nci;              mi = ki/nci;
1254          double          double
1255              *vali = REAL(VECTOR_ELT(ans, i));              *chol = Memcpy(Calloc(ncisq, double),
1256                               REAL(VECTOR_ELT(Omega, i)), ncisq),
1257                *tmp = Calloc(ncisq, double);
1258    
1259          F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);  
1260          if (info)          F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1261              error("DPOTRF returned error code %d for component %d of Omega",          if (j)
1262                    info, i + 1);              error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1263          F77_CALL(dpotri)("U", &nci, vali, &nci, &info);          Memcpy(tmp, chol, ncisq);
1264          if (info)          F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);
1265              error("DPOTRI returned error code %d for component %d of Omega",          if (j)
1266                    info, i + 1);              error("DPOTRI gave error code %d on Omega[[%d]]", j, i + 1);
1267          alpha = (double) -mi;          alpha = (double) -mi;
1268          F77_CALL(dsyrk)("U", "N", &nci, &ki,          F77_CALL(dsyrk)("U", "N", &nci, &ki,
1269                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,
1270                          &alpha, vali, &nci);                          &alpha, tmp, &nci);
1271          alpha = ((double)(REML?(nobs-p):nobs));          alpha = ((double)(REML?(nobs-p):nobs));
1272          F77_CALL(dsyrk)("U", "N", &nci, &mi,          F77_CALL(dsyrk)("U", "N", &nci, &mi,
1273                          &alpha, b + Gp[i], &nci,                          &alpha, b + Gp[i], &nci,
1274                          &one, vali, &nci);                          &one, tmp, &nci);
1275          if (REML) {          if (REML) {
             int j;  
1276              for (j = 0; j < p; j++) {              for (j = 0; j < p; j++) {
1277                  F77_CALL(dsyrk)("U", "N", &nci, &mi,                  F77_CALL(dsyrk)("U", "N", &nci, &mi,
1278                                  &one, RZX + Gp[i] + j*n, &nci,                                  &one, RZX + Gp[i] + j*n, &nci,
1279                                  &one, vali, &nci);                                  &one, tmp, &nci);
1280              }              }
1281          }          }
1282            if (nci == 1) {
1283                REAL(ans)[cind++] = *tmp *
1284                    (uncst ? *REAL(VECTOR_ELT(Omega, i)) : 1.);
1285            } else {
1286                int k, odind = cind + nci;
1287                if (uncst) {
1288                    int ione = 1, kk;
1289                    double *rr = Calloc(nci, double);
1290                    nlme_symmetrize(tmp, nci);
1291                    for (j = 0; j < nci; j++, cind++) {
1292                        for (k = 0; k < nci; k++) rr[k] = chol[j + k*nci];
1293                        REAL(ans)[cind] = 0.;
1294                        for (k = j; k < nci; k++) {
1295                            for (kk = j; kk < nci; kk++) {
1296                                REAL(ans)[cind] += rr[k] * rr[kk] *
1297                                    tmp[kk * nci + k];
1298                            }
1299                        }
1300                        for (k = 0; k < nci; k++) rr[k] *= rr[j];
1301                        for (k = j + 1; k < nci; k++) {
1302                            REAL(ans)[odind++] =
1303                                F77_CALL(ddot)(&nci, rr, &ione, tmp + k, &nci) +
1304                                F77_CALL(ddot)(&nci, rr, &ione,
1305                                               tmp + k*nci, &ione);
1306                        }
1307                    }
1308                    Free(rr);
1309                } else {
1310                    for (j = 0; j < nci; j++) {
1311                        REAL(ans)[cind++] = tmp[j * ncip1];
1312                        for (k = j + 1; k < nci; k++) {
1313                            REAL(ans)[odind++] = tmp[k*nci + j] * 2.;
1314                        }
1315                    }
1316                }
1317                cind = odind;
1318            }
1319            Free(tmp); Free(chol);
1320      }      }
1321      UNPROTECT(1);      UNPROTECT(1);
1322      return ans;      return ans;

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

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