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 111, Tue Apr 20 15:46:04 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 432  Line 432 
432          for (k = j+1; k < nf; k++) { /* off-diagonals */          for (k = j+1; k < nf; k++) { /* off-diagonals */
433              int *fpk = INTEGER(VECTOR_ELT(facs, k)),              int *fpk = INTEGER(VECTOR_ELT(facs, k)),
434                  *Apk = Ap + Gp[k],                  *Apk = Ap + Gp[k],
435                  nck = nc[k];                  nck = nc[k],
436                    scalar = ncj == 1 && nck == 1;
437              double              double
438                  *Zk = Z[k];                  *Zk = Z[k], *work;
439                if (!scalar) work = Calloc(ncj * nck, double);
440              for (i = 0; i < nobs; i++) {              for (i = 0; i < nobs; i++) {
441                  int ii, ind = -1, fpji = fpj[i] - 1,                  int ii, ind = -1, fpji = fpj[i] - 1,
442                      row = Gp[j] + fpji * ncj,                      row = Gp[j] + fpji * ncj,
443                      fpki = fpk[i] - 1,                      fpki = fpk[i] - 1,
444                      lastind = Apk[fpki + 1];                      lastind = Apk[fpki*nck + 1];
445                  for (ii = Apk[fpki]; ii < lastind; ii++) {                  for (ii = Apk[fpki*nck]; ii < lastind; ii++) {
446                      if (Ai[ii] == row) {                      if (Ai[ii] == row) {
447                          ind = ii;                          ind = ii;
448                          break;                          break;
449                      }                      }
450                  }                  }
451                  if (ind < 0) error("logic error in ssclme_update_mm");                  if (ind < 0) error("logic error in ssclme_update_mm");
452                  if (Ncj || nck > 1) {                  if (scalar) {   /* update scalars directly */
453                                  /* FIXME: run a loop to update */                      Ax[ind] += Zj[i] * Zk[i];
454                      error("code not yet written");                  } else {
455                  } else {        /* update scalars directly */                      int jj, offset = ind - Apk[fpki * nck];
456                      Ax[ind] += Zj[fpji] * Zk[fpki];                      F77_CALL(dgemm)("T", "N", &ncj, &nck, &ione, &one,
457                                        Zj + i, &nobs, Zk + i, &nobs,
458                                        &zero, work, &ncj);
459                        for (jj = 0; jj < nck; jj++) {
460                            ind = Apk[fpki * nck + jj] + offset;
461                            if (Ai[ind] != row)
462                                error("logic error in ssclme_update_mm");
463                            for (ii = 0; ii < ncj; ii++) {
464                                Ax[ind++] += work[jj * ncj + ii];
465                  }                  }
466              }              }
467          }          }
468      }      }
469                if (!scalar) Free(work);
470            }
471        }
472      Free(Z);      Free(Z);
473      ssclme_transfer_dimnames(x, facs, mmats);      ssclme_transfer_dimnames(x, facs, mmats);
474      return R_NilValue;      return R_NilValue;
# Line 721  Line 733 
733                          "Rank deficient variance matrix at group %d, level %d",                          "Rank deficient variance matrix at group %d, level %d",
734                          i + 1, j + 1);                          i + 1, j + 1);
735              }              }
736                Free(tmp);
737          }          }
738          return R_NilValue;          return R_NilValue;
739      }      }
# Line 791  Line 804 
804      }      }
805      return R_NilValue;      return R_NilValue;
806  }  }
807    
808  SEXP ssclme_invert(SEXP x)  SEXP ssclme_invert(SEXP x)
809  {  {
810      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
# Line 867  Line 881 
881    
882  /**  /**
883   * Extract the conditional estimates of the fixed effects   * Extract the conditional estimates of the fixed effects
  * FIXME: Add names  
884   *   *
885   * @param x Pointer to an ssclme object   * @param x Pointer to an ssclme object
886   *   *
# Line 895  Line 908 
908    
909  /**  /**
910   * 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.  
911   *   *
912   * @param x Pointer to an ssclme object   * @param x Pointer to an ssclme object
913   *   *
# Line 965  Line 976 
976    
977  /**  /**
978   * Extract the upper triangles of the Omega matrices.   * Extract the upper triangles of the Omega matrices.
979   * (These are not in any sense "coefficients" but the extractor is   * (These aren't "coefficients" but the extractor is
980   * called coef for historical reasons.)   * called coef for historical reasons.)
981   *   *
982   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
# Line 983  Line 994 
994    
995      vind = 0;      vind = 0;
996      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
997          int j, k, nci = nc[i];          int nci = nc[i];
998            if (nci == 1) {
999                vv[vind++] = REAL(VECTOR_ELT(Omega, i))[0];
1000            } else {
1001                int j, k, odind = vind + nci, ncip1 = nci + 1;
1002          double *omgi = REAL(VECTOR_ELT(Omega, i));          double *omgi = REAL(VECTOR_ELT(Omega, i));
1003    
1004          for (j = 0; j < nci; j++) {          for (j = 0; j < nci; j++) {
1005              for (k = 0; k <= j; k++) {                  vv[vind++] = omgi[j * ncip1];
1006                  vv[vind++] = omgi[j*nci + k];                  for (k = j + 1; k < nci; k++) {
1007                        vv[odind++] = omgi[k*nci + j];
1008                    }
1009              }              }
1010                vind = odind;
1011          }          }
1012      }      }
1013      UNPROTECT(1);      UNPROTECT(1);
# Line 996  Line 1015 
1015  }  }
1016    
1017  /**  /**
1018   * Extract the upper triangles of the Omega matrices in the unconstrained   * Extract the unconstrained parameters that determine the
1019   * parameterization.   * Omega matrices. (Called coef for historical reasons.)
  * (These are not in any sense "coefficients" but the extractor is  
  * called coef for historical reasons.)  
1020   *   *
1021   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1022   *   *
1023   * @return numeric vector of the values in the upper triangles of the   * @return numeric vector of unconstrained parameters that determine the
1024   * Omega matrices   * Omega matrices
1025   */   */
1026  SEXP ssclme_coefUnc(SEXP x)  SEXP ssclme_coefUnc(SEXP x)
# Line 1025  Line 1042 
1042                                   REAL(VECTOR_ELT(Omega, i)), ncisq);                                   REAL(VECTOR_ELT(Omega, i)), ncisq);
1043              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
1044              if (j)              /* should never happen */              if (j)              /* should never happen */
1045                  error("DPOTRF returned error code %d", j);                  error("DPOTRF returned error code %d on Omega[[%d]]",
1046                          j, i+1);
1047              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
1048                  double diagj = tmp[j * ncip1];                  double diagj = tmp[j * ncip1];
1049                  vv[vind++] = 2. * log(diagj);                  vv[vind++] = 2. * log(diagj);
# Line 1046  Line 1064 
1064  }  }
1065    
1066  /**  /**
1067   * Assign the upper triangles of the Omega matrices in the unconstrained   * Assign the Omega matrices from the unconstrained parameterization.
  * parameterization.  
1068   *   *
1069   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1070   * @param coef pointer to an numeric vector of appropriate length   * @param coef pointer to an numeric vector of appropriate length
# Line 1078  Line 1095 
1095              double              double
1096                  *omgi = REAL(VECTOR_ELT(Omega, i)),                  *omgi = REAL(VECTOR_ELT(Omega, i)),
1097                  *tmp = Calloc(ncisq, double),                  *tmp = Calloc(ncisq, double),
1098                  diagj, one = 1.;                  diagj, one = 1., zero = 0.;
1099              /* FIXEME: Change this to use a factor and dsyrk */  
                                 /* LD in omgi and L' in tmp */  
1100              memset(omgi, 0, sizeof(double) * ncisq);              memset(omgi, 0, sizeof(double) * ncisq);
1101              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
1102                  omgi[j * ncip1] = diagj = exp(cc[cind++]);                  tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);
1103                  for (k = j + 1; k < nci; k++) {                  for (k = j + 1; k < nci; k++) {
1104                      omgi[j*nci + k] = diagj * (tmp[k*nci + j] = cc[odind++]);                      tmp[k*nci + j] = cc[odind++] * diagj;
1105                  }                  }
1106              }              }
1107              F77_CALL(dtrmm)("R", "U", "N", "U", &nci, &nci, &one,              F77_CALL(dsyrk)("U", "T", &nci, &nci, &one,
1108                              tmp, &nci, omgi, &nci);                              tmp, &nci, &zero, omgi, &nci);
1109              Free(tmp);              Free(tmp);
1110              cind = odind;              cind = odind;
1111          }          }
# Line 1100  Line 1116 
1116    
1117  /**  /**
1118   * Assign the upper triangles of the Omega matrices.   * Assign the upper triangles of the Omega matrices.
1119   * (These are not in any sense "coefficients" but are   * (Called coef for historical reasons.)
  * called coef for historical reasons.)  
1120   *   *
1121   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
1122   * @param coef pointer to an numeric vector of appropriate length   * @param coef pointer to an numeric vector of appropriate length
# Line 1121  Line 1136 
1136                coef_length(nf, nc));                coef_length(nf, nc));
1137      cind = 0;      cind = 0;
1138      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1139          int j, k, nci = nc[i];          int nci = nc[i];
1140            if (nci == 1) {
1141                REAL(VECTOR_ELT(Omega, i))[0] = cc[cind++];
1142            } else {
1143                int j, k, odind = cind + nci, ncip1 = nci + 1;
1144          double *omgi = REAL(VECTOR_ELT(Omega, i));          double *omgi = REAL(VECTOR_ELT(Omega, i));
1145    
1146          for (j = 0; j < nci; j++) {          for (j = 0; j < nci; j++) {
1147              for (k = 0; k <= j; k++) {                  omgi[j * ncip1] = cc[cind++];
1148                  omgi[j*nci + k] = cc[cind++];                  for (k = j + 1; k < nci; k++) {
1149                        omgi[k*nci + j] = cc[odind++];
1150                    }
1151              }              }
1152                cind = odind;
1153          }          }
1154      }      }
1155      status[0] = status[1] = 0;      status[0] = status[1] = 0;
# Line 1209  Line 1232 
1232      return R_NilValue;      return R_NilValue;
1233  }  }
1234    
1235  SEXP ssclme_gradient(SEXP x, SEXP REMLp)  SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)
1236  {  {
1237      SEXP      SEXP
1238            Omega = GET_SLOT(x, Matrix_OmegaSym),
1239          RZXsl = GET_SLOT(x, Matrix_RZXSym),          RZXsl = GET_SLOT(x, Matrix_RZXSym),
         ans = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym))),  
1240          ncsl = GET_SLOT(x, Matrix_ncSym),          ncsl = GET_SLOT(x, Matrix_ncSym),
1241          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
1242      int      int
# Line 1221  Line 1244 
1244          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),          *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1245          *nc = INTEGER(ncsl),          *nc = INTEGER(ncsl),
1246          REML = asLogical(REMLp),          REML = asLogical(REMLp),
1247          i, info,          cind, i, n = dims[0],
1248          n = dims[0],          nf = length(Omega),
1249          nf = length(ncsl) - 2,          nobs, odind, p, pp1 = dims[1],
1250          nobs = nc[nf + 1],          uncst = asLogical(Uncp);
         p,  
         pp1 = dims[1];  
1251      double      double
1252          *RZX = REAL(RZXsl),          *RZX = REAL(RZXsl),
1253          *b,          *b,
1254          alpha,          alpha,
1255          one = 1.;          one = 1.;
1256        SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
1257    
1258        nobs = nc[nf + 1];
1259      p = pp1 - 1;      p = pp1 - 1;
1260      b = RZX + p * n;      b = RZX + p * n;
1261      ssclme_invert(x);      ssclme_invert(x);
1262        cind = 0;
1263      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1264          int ki = Gp[i+1] - Gp[i],          int j, ki = Gp[i+1] - Gp[i],
1265              nci = nc[i],              nci = nc[i], ncip1 = nci + 1, ncisq = nci * nci,
1266              mi = ki/nci;              mi = ki/nci;
1267          double          double
1268              *vali = REAL(VECTOR_ELT(ans, i));              *chol = Memcpy(Calloc(ncisq, double),
1269                               REAL(VECTOR_ELT(Omega, i)), ncisq),
1270                *tmp = Calloc(ncisq, double);
1271    
1272          F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);  
1273          if (info)          F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1274              error("DPOTRF returned error code %d for component %d of Omega",          if (j)
1275                    info, i + 1);              error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1276          F77_CALL(dpotri)("U", &nci, vali, &nci, &info);          Memcpy(tmp, chol, ncisq);
1277          if (info)          F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);
1278              error("DPOTRI returned error code %d for component %d of Omega",          if (j)
1279                    info, i + 1);              error("DPOTRI gave error code %d on Omega[[%d]]", j, i + 1);
1280          alpha = (double) -mi;          alpha = (double) -mi;
1281          F77_CALL(dsyrk)("U", "N", &nci, &ki,          F77_CALL(dsyrk)("U", "N", &nci, &ki,
1282                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,
1283                          &alpha, vali, &nci);                          &alpha, tmp, &nci);
1284          alpha = ((double)(REML?(nobs-p):nobs));          alpha = ((double)(REML?(nobs-p):nobs));
1285          F77_CALL(dsyrk)("U", "N", &nci, &mi,          F77_CALL(dsyrk)("U", "N", &nci, &mi,
1286                          &alpha, b + Gp[i], &nci,                          &alpha, b + Gp[i], &nci,
1287                          &one, vali, &nci);                          &one, tmp, &nci);
1288          if (REML) {          if (REML) {
             int j;  
1289              for (j = 0; j < p; j++) {              for (j = 0; j < p; j++) {
1290                  F77_CALL(dsyrk)("U", "N", &nci, &mi,                  F77_CALL(dsyrk)("U", "N", &nci, &mi,
1291                                  &one, RZX + Gp[i] + j*n, &nci,                                  &one, RZX + Gp[i] + j*n, &nci,
1292                                  &one, vali, &nci);                                  &one, tmp, &nci);
1293                }
1294              }              }
1295            if (nci == 1) {
1296                REAL(ans)[cind++] = *tmp *
1297                    (uncst ? *REAL(VECTOR_ELT(Omega, i)) : 1.);
1298            } else {
1299                int k, odind = cind + nci;
1300                if (uncst) {
1301                    int ione = 1, kk;
1302                    double *rr = Calloc(nci, double);
1303                    nlme_symmetrize(tmp, nci);
1304                    for (j = 0; j < nci; j++, cind++) {
1305                        for (k = 0; k < nci; k++) rr[k] = chol[j + k*nci];
1306                        REAL(ans)[cind] = 0.;
1307                        for (k = j; k < nci; k++) {
1308                            for (kk = j; kk < nci; kk++) {
1309                                REAL(ans)[cind] += rr[k] * rr[kk] *
1310                                    tmp[kk * nci + k];
1311                            }
1312                        }
1313                        for (k = 0; k < nci; k++) rr[k] *= rr[j];
1314                        for (k = j + 1; k < nci; k++) {
1315                            REAL(ans)[odind++] =
1316                                F77_CALL(ddot)(&nci, rr, &ione, tmp + k, &nci) +
1317                                F77_CALL(ddot)(&nci, rr, &ione,
1318                                               tmp + k*nci, &ione);
1319                        }
1320                    }
1321                    Free(rr);
1322                } else {
1323                    for (j = 0; j < nci; j++) {
1324                        REAL(ans)[cind++] = tmp[j * ncip1];
1325                        for (k = j + 1; k < nci; k++) {
1326                            REAL(ans)[odind++] = tmp[k*nci + j] * 2.;
1327                        }
1328                    }
1329                }
1330                cind = odind;
1331          }          }
1332            Free(tmp); Free(chol);
1333      }      }
1334      UNPROTECT(1);      UNPROTECT(1);
1335      return ans;      return ans;

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

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