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

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

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