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 19, Tue Mar 23 23:29:55 2004 UTC revision 21, Sat Mar 27 00:06:02 2004 UTC
# Line 219  Line 219 
219  ssclme_create(SEXP facs, SEXP ncv, SEXP threshold)  ssclme_create(SEXP facs, SEXP ncv, SEXP threshold)
220  {  {
221      SEXP ctab, nms, ssc, tmp,      SEXP ctab, nms, ssc, tmp,
222          val = PROTECT(allocVector(VECSXP, 2));          val = PROTECT(allocVector(VECSXP, 2)),
223            dd = PROTECT(allocVector(INTSXP, 3));   /* dimensions of 3-D arrays */
224      int *Ai, *Ap, *Gp, *LIp, *Lp, *Parent,      int *Ai, *Ap, *Gp, *LIp, *Lp, *Parent,
225          *nc, Lnz, i, nf = length(facs), nzcol, pp1;          *nc, Lnz, i, nf = length(facs), nzcol, pp1,
226            *dims = INTEGER(dd);
227    
228      if (length(ncv) != (nf + 1))      if (length(ncv) != (nf + 1))
229          error("length of nc (%d) should be length of facs (%d) + 1",          error("length of nc (%d) should be length of facs (%d) + 1",
# Line 296  Line 298 
298      tmp = GET_SLOT(ssc, Matrix_bVarSym);      tmp = GET_SLOT(ssc, Matrix_bVarSym);
299      setAttrib(tmp, R_NamesSymbol, getAttrib(facs, R_NamesSymbol));      setAttrib(tmp, R_NamesSymbol, getAttrib(facs, R_NamesSymbol));
300      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
301          int ncol = Gp[i+1] - Gp[i], nrow = nc[i];          int nci = nc[i], mi = (Gp[i+1] - Gp[i])/nc[i];
302    
303          SET_VECTOR_ELT(tmp, i, allocMatrix(REALSXP, nrow, ncol));          dims[0] = dims[1] = nci;
304            dims[2] = mi;
305            SET_VECTOR_ELT(tmp, i, allocArray(REALSXP, dd));
306          memset(REAL(VECTOR_ELT(tmp, i)), 0,          memset(REAL(VECTOR_ELT(tmp, i)), 0,
307                 sizeof(double) * nrow * ncol);                 sizeof(double) * nci * nci * mi);
308      }      }
309      SET_SLOT(ssc, Matrix_LIpSym, allocVector(INTSXP, nzcol + 1));      SET_SLOT(ssc, Matrix_LIpSym, allocVector(INTSXP, nzcol + 1));
310      LIp = INTEGER(GET_SLOT(ssc, Matrix_LIpSym));      LIp = INTEGER(GET_SLOT(ssc, Matrix_LIpSym));
# Line 313  Line 317 
317          memset(REAL(GET_SLOT(ssc, Matrix_LIxSym)), 0,          memset(REAL(GET_SLOT(ssc, Matrix_LIxSym)), 0,
318                 sizeof(double) * Lnz);                 sizeof(double) * Lnz);
319      }      }
320      UNPROTECT(1);      UNPROTECT(2);
321      return val;      return val;
322  }  }
323    
# Line 876  Line 880 
880   */   */
881  SEXP ssclme_ranef(SEXP x)  SEXP ssclme_ranef(SEXP x)
882  {  {
883      SEXP RZXsl = GET_SLOT(x, Matrix_RZXSym);      SEXP RZXsl = GET_SLOT(x, Matrix_RZXSym),
884            GpSl = GET_SLOT(x, Matrix_GpSym);
885      int *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),      int *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
886          j,          *Gp = INTEGER(GpSl),
887            *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
888            i, j,
889          n = dims[0],          n = dims[0],
890            nf = length(GpSl) - 1,
891          pp1 = dims[1],          pp1 = dims[1],
892          p = pp1 - 1;          p = pp1 - 1;
893      SEXP val = PROTECT(allocVector(REALSXP, n));      SEXP val = PROTECT(allocVector(VECSXP, nf));
894      double      double
895          *RZX = REAL(RZXsl),          *b = REAL(RZXsl) + n * p,
         *b = REAL(val),  
896          ryyinv;         /* ryy-inverse */          ryyinv;         /* ryy-inverse */
897    
898      ssclme_invert(x);      ssclme_invert(x);
     Memcpy(b, RZX + p * n, n);  
899      ryyinv = REAL(GET_SLOT(x, Matrix_RXXSym))[pp1*pp1 - 1];      ryyinv = REAL(GET_SLOT(x, Matrix_RXXSym))[pp1*pp1 - 1];
900      for (j = 0; j < n; j++) b[j] /= ryyinv;      for (i = 0; i < nf; i++) {
901            int nci = nc[i], Mi = (Gp[i+1] - Gp[i]), mi = Mi/nci;
902            double *mm;
903    
904            SET_VECTOR_ELT(val, i, allocMatrix(REALSXP, nci, mi));
905            mm = Memcpy(REAL(VECTOR_ELT(val, i)), b, Mi);
906            b += nci * mi;
907            for (j = 0; j < Mi; j++) mm[j] /= ryyinv;
908        }
909      UNPROTECT(1);      UNPROTECT(1);
910      return val;      return val;
911  }  }

Legend:
Removed from v.19  
changed lines
  Added in v.21

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