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 27, Sun Mar 28 16:40:05 2004 UTC revision 28, Sun Mar 28 17:06:05 2004 UTC
# Line 1080  Line 1080 
1080                  *omgi = REAL(VECTOR_ELT(Omega, i)),                  *omgi = REAL(VECTOR_ELT(Omega, i)),
1081                  *tmp = Calloc(ncisq, double),                  *tmp = Calloc(ncisq, double),
1082                  diagj, one = 1.;                  diagj, one = 1.;
1083                /* FIXEME: Change this to use a factor and dsyrk */
1084                                  /* LD in omgi and L' in tmp */                                  /* LD in omgi and L' in tmp */
1085              memset(omgi, 0, sizeof(double) * ncisq);              memset(omgi, 0, sizeof(double) * ncisq);
1086              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
# Line 1236  Line 1237 
1237      return val;      return val;
1238  }  }
1239    
1240    
1241    SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats)
1242    {
1243        SEXP val, b;
1244        int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1245            i, ione = 1, nf = length(facs), nobs, p;
1246        double *vv, one = 1.0, zero = 0.0;
1247    
1248        if (nf < 1)
1249            error("null factor list passed to ssclme_fitted");
1250        nobs = length(VECTOR_ELT(facs, 0));
1251        val = PROTECT(allocVector(REALSXP, nobs));
1252        vv = REAL(val);
1253        p = nc[nf] - 1;
1254        if (p > 0) {
1255            F77_CALL(dgemm)("N", "N", &nobs, &ione, &p, &one,
1256                            REAL(VECTOR_ELT(mmats, nf)), &nobs,
1257                            REAL(PROTECT(ssclme_fixef(x))), &p,
1258                            &zero, vv, &nobs);
1259            UNPROTECT(1);
1260        } else {
1261            memset(vv, 0, sizeof(double) * nobs);
1262        }
1263        b = PROTECT(ssclme_ranef(x));
1264        for (i = 0; i < nf; i++) {
1265            int *ff = INTEGER(VECTOR_ELT(facs, i)), j, nci = nc[i];
1266            double *bb = REAL(VECTOR_ELT(b, i)),
1267                *mm = REAL(VECTOR_ELT(mmats, i));
1268            for (j = 0; j < nobs; ) {
1269                int nn = 1, lev = ff[j];
1270                /* check for adjacent rows with same factor level */
1271                while (ff[j + nn] == lev) nn++;
1272                F77_CALL(dgemm)("N", "N", &nn, &ione, &nci,
1273                                &one, mm + j, &nobs,
1274                                bb + lev - 1, &nci,
1275                                &one, vv + j, &nobs);
1276                j += nn;
1277            }
1278        }
1279        UNPROTECT(2);
1280        return val;
1281    }

Legend:
Removed from v.27  
changed lines
  Added in v.28

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