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 115, Wed Apr 21 22:37:28 2004 UTC revision 125, Sun Apr 25 11:59:37 2004 UTC
# Line 601  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("Could not factor downdated X'X, 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 612  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 809  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 957  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 1258  Line 1263 
1263          alpha,          alpha,
1264          one = 1.;          one = 1.;
1265      SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));      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;
# Line 1381  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 1406  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;
     UNPROTECT(2);  
     return val;  
1426  }  }
1427    
1428  SEXP ssclme_collapse(SEXP x)  SEXP ssclme_collapse(SEXP x)
# Line 1437  Line 1445 
1445    
1446      INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;      INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;
1447      SET_SLOT(ans, Matrix_XtXSym, allocMatrix(REALSXP, 1, 1));      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));      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));      SET_SLOT(ans, Matrix_ZtXSym, allocMatrix(REALSXP, nz, 1));
1452      SET_SLOT(ans, Matrix_RZXSym, allocMatrix(REALSXP, nz, 1));      SET_SLOT(ans, Matrix_RZXSym, allocMatrix(REALSXP, nz, 1));
1453      LOGICAL(GET_SLOT(ans, Matrix_statusSym))[0] = 0;      LOGICAL(GET_SLOT(ans, Matrix_statusSym))[0] = 0;

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