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 105, Sun Apr 18 22:12:10 2004 UTC
# Line 1224  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),
         ans,  
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 1233  Line 1232 
1232          *nc = INTEGER(ncsl),          *nc = INTEGER(ncsl),
1233          REML = asLogical(REMLp),          REML = asLogical(REMLp),
1234          cind, i, n = dims[0],          cind, i, n = dims[0],
1235          nf = length(ncsl) - 2,          nf = length(Omega),
1236          nobs, odind, p, pp1 = dims[1],          nobs, odind, p, pp1 = dims[1],
1237          uncst = asLogical(Uncp);          uncst = asLogical(Uncp);
1238      double      double
# Line 1241  Line 1240 
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];      nobs = nc[nf + 1];
1246      p = pp1 - 1;      p = pp1 - 1;
1247      b = RZX + p * n;      b = RZX + p * n;
     ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));  
1248      ssclme_invert(x);      ssclme_invert(x);
1249      cind = 0;      cind = 0;
1250      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1251          int j, 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              *tmp = Memcpy(Calloc(nci * nci, double),              *chol = Memcpy(Calloc(ncisq, double),
1256                            REAL(VECTOR_ELT(Omega, i)), nci * nci);                             REAL(VECTOR_ELT(Omega, i)), ncisq),
1257                *tmp = Calloc(ncisq, double);
1258    
1259          F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);  
1260            F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1261          if (j)          if (j)
1262              error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);              error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
1263            Memcpy(tmp, chol, ncisq);
1264          F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);          F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);
1265          if (j)          if (j)
1266              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 1283 
1283              REAL(ans)[cind++] = *tmp *              REAL(ans)[cind++] = *tmp *
1284                  (uncst ? *REAL(VECTOR_ELT(Omega, i)) : 1.);                  (uncst ? *REAL(VECTOR_ELT(Omega, i)) : 1.);
1285          } else {          } else {
1286                int odind = cind + nci;
1287              if (uncst) {              if (uncst) {
1288                  error("Code not written yet");                  int kk;
1289                    nlme_symmetrize(tmp, nci);
1290                    for (j = 0; j < nci; j++, cind++) {
1291                        double *dder = REAL(ans) + cind;
1292                        *dder = 0.;
1293                        for (k = j; k < nci; k++) {
1294                            for (kk = j; kk < nci; kk++) {
1295                                *dder += chol[k*nci + j] * chol[kk*nci + j] *
1296                                    tmp[kk * nci + k];
1297                            }
1298                        }
1299                        for (k = j + 1; k < nci; k++) {
1300                            REAL(ans)[odind++] = -1.;
1301                        }
1302                    }
1303              } else {              } else {
                 int k, ncip1 = nci + 1, odind = cind + nci;  
1304                  for (j = 0; j < nci; j++) {                  for (j = 0; j < nci; j++) {
1305                      REAL(ans)[cind++] = tmp[j * ncip1];                      REAL(ans)[cind++] = tmp[j * ncip1];
1306                      for (k = j + 1; k < nci; k++) {                      for (k = j + 1; k < nci; k++) {
1307                          REAL(ans)[odind++] = tmp[k*nci + j] * 2.;                          REAL(ans)[odind++] = tmp[k*nci + j] * 2.;
1308                      }                      }
1309                  }                  }
                 cind = odind;  
1310              }              }
1311                cind = odind;
1312          }          }
1313          Free(tmp);          Free(tmp); Free(chol);
1314      }      }
1315      UNPROTECT(1);      UNPROTECT(1);
1316      return ans;      return ans;

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

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