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 104, Sun Apr 18 19:21:57 2004 UTC revision 125, Sun Apr 25 11:59:37 2004 UTC
# 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);
# Line 589  Line 601 
601          F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,          F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,
602                          RZX, &n, &one, RXX, &pp1);                          RZX, &n, &one, RXX, &pp1);
603          F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);          F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);
604          if (i)          if (i) {
605              error("DPOTRF returned error code %d", i);              warning("Could not factor downdated X'X, code %d", i);
606                dcmp[2] = dcmp[3] = deviance[0] = deviance[1] = NA_REAL;
607            } else {
608                                  /* logdet of RXX */                                  /* logdet of RXX */
609          for (i = 0; i < (pp1 - 1); i++)          for (i = 0; i < (pp1 - 1); i++)
610              dcmp[2] += 2 * log(RXX[i*pp2]);              dcmp[2] += 2 * log(RXX[i*pp2]);
# Line 600  Line 614 
614              dcmp[0] - dcmp[1] + nobs*(1+dcmp[3]+log(2*PI/nobs));              dcmp[0] - dcmp[1] + nobs*(1+dcmp[3]+log(2*PI/nobs));
615          deviance[1] = dcmp[0] - dcmp[1] + /* REML */          deviance[1] = dcmp[0] - dcmp[1] + /* REML */
616              dcmp[2] + nreml * (1. + dcmp[3] + log(2. * PI/nreml));              dcmp[2] + nreml * (1. + dcmp[3] + log(2. * PI/nreml));
617            }
618          status[0] = 1;          status[0] = 1;
619          status[1] = 0;          status[1] = 0;
620      }      }
# Line 696  Line 711 
711              }              }
712              tmp = Calloc(nr * nci, double); /* scratch storage */              tmp = Calloc(nr * nci, double); /* scratch storage */
713              nr1 = nr + 1;              nr1 = nr + 1;
714                                  /* initialize bVi to zero (cosmetic) */                                  /* initialize bVi to zero */
715              memset(bVi, 0, sizeof(double) * (G2 - G1) * nci);              memset(bVi, 0, sizeof(double) * (G2 - G1) * nci);
716              for (j = G1; j < G2; j += nci) {              for (j = G1; j < G2; j += nci) {
717                  memset(tmp, 0, sizeof(double) * nr * nci);                  memset(tmp, 0, sizeof(double) * nr * nci);
# Line 712  Line 727 
727                          tmp[k * nr + kk] *= DIsqrt[LIi[LIp[j] + kk - 1]];                          tmp[k * nr + kk] *= DIsqrt[LIi[LIp[j] + kk - 1]];
728                      }                      }
729                  }                  }
730                  F77_CALL(dsyrk)("U", "T", &nci, &rr, &one, tmp, &nr,                  F77_CALL(dsyrk)("L", "T", &nci, &rr, &one, tmp, &nr,
731                                  &zero, bVi + (j - G1) * nci, &nci);                                  &zero, bVi + (j - G1) * nci, &nci);
732                  F77_CALL(dpotrf)("U", &nci, bVi + (j - G1) * nci,                  F77_CALL(dpotrf)("L", &nci, bVi + (j - G1) * nci,
733                                   &nci, &kk);                                   &nci, &kk);
734                  if (kk)         /* should never happen */                  if (kk)         /* should never happen */
735                      error(                      error(
736                          "Rank deficient variance matrix at group %d, level %d",                          "Rank deficient variance matrix at group %d, level %d",
737                          i + 1, j + 1);                          i + 1, j + 1);
738              }              }
739                Free(tmp);
740          }          }
741          return R_NilValue;          return R_NilValue;
742      }      }
# Line 777  Line 793 
793                          tmp[k * nr + kk] *= DIsqrt[ind[kk]];                          tmp[k * nr + kk] *= DIsqrt[ind[kk]];
794                      }                      }
795                  }                  }
796                  F77_CALL(dsyrk)("U", "T", &nci, &nr, &one, tmp, &nr,                  F77_CALL(dsyrk)("L", "T", &nci, &nr, &one, tmp, &nr,
797                                  &zero, mpt + (j - Gp[i])*nci, &nci);                                  &zero, mpt + (j - Gp[i])*nci, &nci);
798                  F77_CALL(dpotrf)("U", &nci, mpt + (j - Gp[i])*nci,                  F77_CALL(dpotrf)("L", &nci, mpt + (j - Gp[i])*nci,
799                                   &nci, &info);                                   &nci, &info);
800                  if (info)       /* should never happen */                  if (info)       /* should never happen */
801                      error(                      error(
# Line 796  Line 812 
812  {  {
813      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
814      if (!status[0]) ssclme_factor(x);      if (!status[0]) ssclme_factor(x);
815        if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0]))
816            error("Unable to invert singular factor of downdated X'X");
817      if (!status[1]) {      if (!status[1]) {
818          SEXP          SEXP
819              RZXsl = GET_SLOT(x, Matrix_RZXSym);              RZXsl = GET_SLOT(x, Matrix_RZXSym);
# Line 944  Line 962 
962  {  {
963      SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);      SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);
964      int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1],      int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1],
965          nobs = INTEGER(GET_SLOT(x, Matrix_ncSym))          nobs = INTEGER(GET_SLOT(x, Matrix_ncSym))[
966          [length(GET_SLOT(x, Matrix_GpSym))];              length(GET_SLOT(x, Matrix_GpSym))];
967    
968      ssclme_invert(x);      ssclme_invert(x);
969      return ScalarReal(1./(REAL(RXXsl)[pp1*pp1 - 1] *      return ScalarReal(1./(REAL(RXXsl)[pp1*pp1 - 1] *
# Line 1209  Line 1227 
1227                  }                  }
1228              }              }
1229              F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);              F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);
1230              if (info) error("DPOTRF returned error code %d", info);              if (info)
1231                    error("DPOTRF returned error code %d in Omega[[%d]] update",
1232                          info, i + 1);
1233              F77_CALL(dpotri)("U", &nci, vali, &nci, &info);              F77_CALL(dpotri)("U", &nci, vali, &nci, &info);
1234              if (info) error("DPOTRI returned error code %d", info);              if (info)
1235                    error("DPOTRI returned error code %d in Omega[[%d]] update",
1236                          info, i + 1);
1237          }          }
1238          status[0] = status[1] = 0;          status[0] = status[1] = 0;
1239      }      }
# Line 1224  Line 1246 
1246      SEXP      SEXP
1247          Omega = GET_SLOT(x, Matrix_OmegaSym),          Omega = GET_SLOT(x, Matrix_OmegaSym),
1248          RZXsl = GET_SLOT(x, Matrix_RZXSym),          RZXsl = GET_SLOT(x, Matrix_RZXSym),
         ans,  
1249          ncsl = GET_SLOT(x, Matrix_ncSym),          ncsl = GET_SLOT(x, Matrix_ncSym),
1250          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
1251      int      int
# Line 1233  Line 1254 
1254          *nc = INTEGER(ncsl),          *nc = INTEGER(ncsl),
1255          REML = asLogical(REMLp),          REML = asLogical(REMLp),
1256          cind, i, n = dims[0],          cind, i, n = dims[0],
1257          nf = length(ncsl) - 2,          nf = length(Omega),
1258          nobs, odind, p, pp1 = dims[1],          nobs, p, pp1 = dims[1],
1259          uncst = asLogical(Uncp);          uncst = asLogical(Uncp);
1260      double      double
1261          *RZX = REAL(RZXsl),          *RZX = REAL(RZXsl),
1262          *b,          *b,
1263          alpha,          alpha,
1264          one = 1.;          one = 1.;
1265        SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
1266        int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
1267    
1268        ssclme_factor(x);
1269        if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0])) {
1270            int ncoef = coef_length(nf, nc);
1271            for (i = 0; i < ncoef; i++) REAL(ans)[i] = NA_REAL;
1272            UNPROTECT(1);
1273            return ans;
1274        }
1275      nobs = nc[nf + 1];      nobs = nc[nf + 1];
1276      p = pp1 - 1;      p = pp1 - 1;
1277      b = RZX + p * n;      b = RZX + p * n;
     ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));  
1278      ssclme_invert(x);      ssclme_invert(x);
1279      cind = 0;      cind = 0;
1280      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1281          int j, ki = Gp[i+1] - Gp[i],          int j, ki = Gp[i+1] - Gp[i],
1282              nci = nc[i],              nci = nc[i], ncip1 = nci + 1, ncisq = nci * nci,
1283              mi = ki/nci;              mi = ki/nci;
1284          double          double
1285              *tmp = Memcpy(Calloc(nci * nci, double),              *chol = Memcpy(Calloc(ncisq, double),
1286                            REAL(VECTOR_ELT(Omega, i)), nci * nci);                             REAL(VECTOR_ELT(Omega, i)), ncisq),
1287                *tmp = Calloc(ncisq, double);
1288    
1289          F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);  
1290            F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1291          if (j)          if (j)
1292              error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);              error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1293            Memcpy(tmp, chol, ncisq);
1294          F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);          F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);
1295          if (j)          if (j)
1296              error("DPOTRI gave error code %d on Omega[[%d]]", j, i + 1);              error("DPOTRI gave error code %d on Omega[[%d]]", j, i + 1);
# Line 1281  Line 1313 
1313              REAL(ans)[cind++] = *tmp *              REAL(ans)[cind++] = *tmp *
1314                  (uncst ? *REAL(VECTOR_ELT(Omega, i)) : 1.);                  (uncst ? *REAL(VECTOR_ELT(Omega, i)) : 1.);
1315          } else {          } else {
1316                int k, odind = cind + nci;
1317              if (uncst) {              if (uncst) {
1318                  error("Code not written yet");                  int ione = 1, kk;
1319                    double *rr = Calloc(nci, double);
1320                    nlme_symmetrize(tmp, nci);
1321                    for (j = 0; j < nci; j++, cind++) {
1322                        for (k = 0; k < nci; k++) rr[k] = chol[j + k*nci];
1323                        REAL(ans)[cind] = 0.;
1324                        for (k = j; k < nci; k++) {
1325                            for (kk = j; kk < nci; kk++) {
1326                                REAL(ans)[cind] += rr[k] * rr[kk] *
1327                                    tmp[kk * nci + k];
1328                            }
1329                        }
1330                        for (k = 0; k < nci; k++) rr[k] *= rr[j];
1331                        for (k = j + 1; k < nci; k++) {
1332                            REAL(ans)[odind++] =
1333                                F77_CALL(ddot)(&nci, rr, &ione, tmp + k, &nci) +
1334                                F77_CALL(ddot)(&nci, rr, &ione,
1335                                               tmp + k*nci, &ione);
1336                        }
1337                    }
1338                    Free(rr);
1339              } else {              } else {
                 int k, ncip1 = nci + 1, odind = cind + nci;  
1340                  for (j = 0; j < nci; j++) {                  for (j = 0; j < nci; j++) {
1341                      REAL(ans)[cind++] = tmp[j * ncip1];                      REAL(ans)[cind++] = tmp[j * ncip1];
1342                      for (k = j + 1; k < nci; k++) {                      for (k = j + 1; k < nci; k++) {
1343                          REAL(ans)[odind++] = tmp[k*nci + j] * 2.;                          REAL(ans)[odind++] = tmp[k*nci + j] * 2.;
1344                      }                      }
1345                  }                  }
                 cind = odind;  
1346              }              }
1347                cind = odind;
1348          }          }
1349          Free(tmp);          Free(tmp); Free(chol);
1350      }      }
1351      UNPROTECT(1);      UNPROTECT(1);
1352      return ans;      return ans;
# Line 1342  Line 1394 
1394      return val;      return val;
1395  }  }
1396    
1397  SEXP ssclme_variances(SEXP x, SEXP REML)  /**
1398     * Return the unscaled variances
1399     *
1400     * @param x pointer to an ssclme object
1401     *
1402     * @return a list similar to the Omega list with the unscaled variances
1403     */
1404    SEXP ssclme_variances(SEXP x)
1405  {  {
1406      SEXP Omg = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym))),      SEXP Omg = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym)));
         val = PROTECT(allocVector(VECSXP, 2));  
1407      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1408          i, nf = length(Omg);          i, nf = length(Omg);
     double sigmasq;  
1409    
   
     SET_VECTOR_ELT(val, 0, Omg);  
     SET_VECTOR_ELT(val, 1, ssclme_sigma(x, REML));  
     sigmasq = REAL(VECTOR_ELT(val, 1))[0];  
     sigmasq = (sigmasq) * (sigmasq);  
1410      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1411          double *mm = REAL(VECTOR_ELT(Omg, i));          double *mm = REAL(VECTOR_ELT(Omg, i));
1412          int j, k, nci = nc[i], ncip1 = nci+1;          int j, nci = nc[i];
1413    
1414          F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);          F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);
1415          if (j)                  /* shouldn't happen */          if (j)                  /* shouldn't happen */
# Line 1367  Line 1419 
1419          if (j)                  /* shouldn't happen */          if (j)                  /* shouldn't happen */
1420              error("DTRTRI returned error code %d on Omega[%d]",              error("DTRTRI returned error code %d on Omega[%d]",
1421                    j, i + 1);                    j, i + 1);
1422          for (j = 0; j < nci; j++) {          nlme_symmetrize(mm, nci);
             mm[j * ncip1] *= sigmasq;  
             for (k = 0; k < j; k++) {  
                 mm[j + k * nci] = (mm[k + j * nci] *= sigmasq);  
             }  
1423          }          }
1424        UNPROTECT(1);
1425        return Omg;
1426      }      }
1427      UNPROTECT(2);  
1428      return val;  SEXP ssclme_collapse(SEXP x)
1429    {
1430        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("ssclme"))),
1431            Omega = GET_SLOT(x, Matrix_OmegaSym),
1432            Dim = GET_SLOT(x, Matrix_DimSym);
1433        int i, nf = length(Omega), nz = INTEGER(Dim)[1];
1434        SEXP copy[] = {Matrix_DSym, Matrix_DIsqrtSym, Matrix_DimSym,
1435                       Matrix_GpSym, Matrix_LIiSym, Matrix_LIpSym,
1436                       Matrix_LIxSym, Matrix_LiSym, Matrix_LpSym,
1437                       Matrix_LxSym, Matrix_OmegaSym, Matrix_ParentSym,
1438                       Matrix_LIxSym, Matrix_LiSym, Matrix_LpSym,
1439                       Matrix_bVarSym, Matrix_devianceSym,
1440                       Matrix_devCompSym, Matrix_iSym, Matrix_ncSym,
1441                       Matrix_statusSym, Matrix_pSym, Matrix_xSym};
1442    
1443        for (i = 0; i < 23; i++)
1444            SET_SLOT(ans, copy[i], duplicate(GET_SLOT(x, copy[i])));
1445    
1446        INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;
1447        SET_SLOT(ans, Matrix_XtXSym, allocMatrix(REALSXP, 1, 1));
1448        REAL(GET_SLOT(ans, Matrix_XtXSym))[0] = NA_REAL;
1449        SET_SLOT(ans, Matrix_RXXSym, allocMatrix(REALSXP, 1, 1));
1450        REAL(GET_SLOT(ans, Matrix_RXXSym))[0] = NA_REAL;
1451        SET_SLOT(ans, Matrix_ZtXSym, allocMatrix(REALSXP, nz, 1));
1452        SET_SLOT(ans, Matrix_RZXSym, allocMatrix(REALSXP, nz, 1));
1453        LOGICAL(GET_SLOT(ans, Matrix_statusSym))[0] = 0;
1454        UNPROTECT(1);
1455        return ans;
1456  }  }
1457    
1458    

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

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