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 45, Mon Apr 5 01:06:55 2004 UTC revision 57, Tue Apr 6 22:09:52 2004 UTC
# Line 1253  Line 1253 
1253    
1254  SEXP ssclme_variances(SEXP x, SEXP REML)  SEXP ssclme_variances(SEXP x, SEXP REML)
1255  {  {
1256      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym),      SEXP Omg = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym))),
1257          val = PROTECT(allocVector(VECSXP, length(Omega) + 1));          val = PROTECT(allocVector(VECSXP, 2));
1258      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),      int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1259          i, nf = length(Omega);          i, nf = length(Omg);
1260      double *sig, sigmasq;      double sigmasq;
1261    
1262      SET_VECTOR_ELT(val, nf, ssclme_sigma(x, REML));  
1263      sig = REAL(VECTOR_ELT(val, nf));      SET_VECTOR_ELT(val, 0, Omg);
1264      sigmasq = (*sig) * (*sig);      SET_VECTOR_ELT(val, 1, ssclme_sigma(x, REML));
1265      *sig = sigmasq;      sigmasq = REAL(VECTOR_ELT(val, 1))[0];
1266        sigmasq = (sigmasq) * (sigmasq);
1267      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1268          SEXP vari = duplicate(VECTOR_ELT(Omega, i));          double *mm = REAL(VECTOR_ELT(Omg, i));
         double *mm = REAL(vari);  
1269          int j, k, nci = nc[i], ncip1 = nci+1;          int j, k, nci = nc[i], ncip1 = nci+1;
1270    
         SET_VECTOR_ELT(val, i, vari);  
1271          F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);          F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);
1272          if (j)                  /* shouldn't happen */          if (j)                  /* shouldn't happen */
1273              error("DPOTRF returned error code %d on Omega[%d]",              error("DPOTRF returned error code %d on Omega[%d]",
# Line 1284  Line 1283 
1283              }              }
1284          }          }
1285      }      }
1286      UNPROTECT(1);      UNPROTECT(2);
1287      return val;      return val;
1288  }  }

Legend:
Removed from v.45  
changed lines
  Added in v.57

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