SCM

SCM Repository

[matrix] Diff of /pkg/src/lmer.c
ViewVC logotype

Diff of /pkg/src/lmer.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 421, Mon Jan 3 20:10:42 2005 UTC revision 422, Tue Jan 4 23:45:15 2005 UTC
# Line 513  Line 513 
513                  int ntot = xdims[0] * xdims[1];                  int ntot = xdims[0] * xdims[1];
514    
515                  /* ZZpO_{j,i} := ZZpO_{j,i} %*% Linv_{i,i}^T %*% D_i^{-1/2} */                  /* ZZpO_{j,i} := ZZpO_{j,i} %*% Linv_{i,i}^T %*% D_i^{-1/2} */
516                  /* FIXME: This should be in L, not ZZpO or change L to diagonal blocks*/                  /* FIXME: Change L to diagonal blocks - off-diagonals are not used*/
517                  cscb_trcbm('R', 'L', 'T', 'U', 1.0, LinvP, ZZOji);                  cscb_trcbm('R', 'L', 'T', 'U', 1.0, LinvP, ZZOji);
518                  for (jj = 0; jj < nlev; jj++) {                  for (jj = 0; jj < nlev; jj++) {
519                      int k, k2 = ZZp[jj + 1];                      int k, k2 = ZZp[jj + 1];
520                      for (k = ZZp[jj]; k < k2; k++)                      for (k = ZZp[jj]; k < k2; k++)
521                          F77_CALL(dtrsm)("R", "U", "N", "N", &nci, xdims + 1,                          F77_CALL(dtrsm)("R", "U", "N", "N", xdims, xdims + 1,
522                                          &one, D + jj * ncisqr, &nci,                                          &one, D + jj * ncisqr, &nci,
523                                          ZZO + k * ntot, xdims);                                          ZZO + k * ntot, xdims);
524                  }                  }
# Line 548  Line 548 
548                  for (jj = 0; jj < nlev; jj++) {                  for (jj = 0; jj < nlev; jj++) {
549                      int k, k2 = ZZp[jj + 1];                      int k, k2 = ZZp[jj + 1];
550                      for (k = ZZp[jj]; k < k2; k++)                      for (k = ZZp[jj]; k < k2; k++)
551                          F77_CALL(dtrsm)("R", "U", "T", "N", &nci, xdims + 1,                          F77_CALL(dtrsm)("R", "U", "T", "N", xdims, xdims + 1,
552                                          &one, D + jj * ncisqr, &nci,                                          &one, D + jj * ncisqr, &nci,
553                                          ZZO + k * ntot, xdims);                                          ZZO + k * ntot, xdims);
554                  }                  }
# Line 641  Line 641 
641              *RZX = REAL(RZXsl),              *RZX = REAL(RZXsl),
642               minus1 = -1., one = 1.;               minus1 = -1., one = 1.;
643    
644            /* RXX := RXX^{-1} */
645          F77_CALL(dtrtri)("U", "N", dims + 1, RXX, dims + 1, &i);          F77_CALL(dtrtri)("U", "N", dims + 1, RXX, dims + 1, &i);
646          if (i)          if (i)
647              error("Leading minor of size %d of downdated X'X,is indefinite",              error("Leading minor of size %d of downdated X'X,is indefinite",
648                      i + 1);                      i + 1);
649            /* RZX := - RZX %*% RXX */
650          F77_CALL(dtrmm)("R", "U", "N", "N", dims, dims + 1, &minus1,          F77_CALL(dtrmm)("R", "U", "N", "N", dims, dims + 1, &minus1,
651                          RXX, dims + 1, RZX, dims);                          RXX, dims + 1, RZX, dims);
652          for(i = 0; i < nf; i++) {          for(i = 0; i < nf; i++) {
# Line 653  Line 655 
655              double *Di = REAL(VECTOR_ELT(Dsl, i)),              double *Di = REAL(VECTOR_ELT(Dsl, i)),
656                  *RZXi = RZX + Gp[i];                  *RZXi = RZX + Gp[i];
657    
658                /* D_i := D_i^{-1}; RZX_i := D_i %*% RZX_i */
659              if (nci == 1) {              if (nci == 1) {
660                  for (j = 0; j < nlev; j++) {                  for (j = 0; j < nlev; j++) {
661                      Di[j] = 1./Di[j];                      Di[j] = 1./Di[j];
# Line 669  Line 672 
672                  }                  }
673              }              }
674          }          }
675          lmer_sm('L', 'T', nf, Gp, dims[1], -1.0, GET_SLOT(x, Matrix_ZZpOSym),          /* RZX := L^{-T} %*% RZX */
676            lmer_sm('L', 'T', nf, Gp, dims[1], 1.0, GET_SLOT(x, Matrix_ZZpOSym),
677                  GET_SLOT(x, Matrix_LinvSym), RZX, dims[0]);                  GET_SLOT(x, Matrix_LinvSym), RZX, dims[0]);
678          status[1] = 1;          status[1] = 1;
679      }      }
# Line 883  Line 887 
887      double      double
888          *b = REAL(RZXsl) + dims[0] * (dims[1] - 1),          *b = REAL(RZXsl) + dims[0] * (dims[1] - 1),
889          nryyinv;                /* negative ryy-inverse */          nryyinv;                /* negative ryy-inverse */
     /* FIXME: Check if this really should be negative ryy-inverse now  
      * that the negative sign is used in the inversion of the RZX slot */  
890    
891      lmer_invert(x);      lmer_invert(x);
892      setAttrib(val, R_NamesSymbol,      setAttrib(val, R_NamesSymbol,
# Line 926  Line 928 
928  static  static
929  SEXP lmer_firstDer(SEXP x, SEXP val)  SEXP lmer_firstDer(SEXP x, SEXP val)
930  {  {
931      SEXP Omega = GET_SLOT(x, Matrix_OmegaSym),      SEXP D = PROTECT(duplicate(GET_SLOT(x, Matrix_DSym))),
932          D = GET_SLOT(x, Matrix_DSym),          LinvP = GET_SLOT(x, Matrix_LinvSym),
933            ZZOP = GET_SLOT(x, Matrix_ZZpOSym),
934            Omega = GET_SLOT(x, Matrix_OmegaSym),
935          RZXsl = GET_SLOT(x, Matrix_RZXSym);          RZXsl = GET_SLOT(x, Matrix_RZXSym);
936      int *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),      int *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
937          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
# Line 936  Line 940 
940          *b = REAL(RZXsl) + dims[0] * p;          *b = REAL(RZXsl) + dims[0] * p;
941    
942      lmer_invert(x);      lmer_invert(x);
943      setAttrib(val, R_NamesSymbol,      for (i = nf - 1; i >= 0; i++) {
               duplicate(getAttrib(Omega, R_NamesSymbol)));  
     for (i = 0; i < nf; i++) {  
944          SEXP DiP = VECTOR_ELT(D, i);          SEXP DiP = VECTOR_ELT(D, i);
945          int *ddims = INTEGER(getAttrib(DiP, R_DimSymbol)), j;          int *ddims = INTEGER(getAttrib(DiP, R_DimSymbol)), j;
946          int nci = ddims[0], nlev = ddims[2];          int nci = ddims[0], nlev = ddims[2];
947          int ncisqr = nci * nci, RZXrows = Gp[i + 1] - Gp[i];          int ncisqr = nci * nci, RZXrows = Gp[i + 1] - Gp[i];
948          double *Di = REAL(DiP), *RZXi = RZX + Gp[i],          double *Di = REAL(DiP), *RZXi = RZX + Gp[i],
949              *bi = b + Gp[i], *mm, dlev = (double) nlev;              *bi = b + Gp[i], *mm = REAL(VECTOR_ELT(val, i)),
950                *tmp = Memcpy(Calloc(ncisqr, double),
951                              REAL(VECTOR_ELT(Omega, i)), ncisqr),
952                dlev = (double) nlev,
953                one = 1., zero = 0.;
954    
         mm = memset(REAL(VECTOR_ELT(val, i)), 0, sizeof(double) * ncisqr * 4);  
955          if (nci == 1) {          if (nci == 1) {
956              int ione = 1;              int ione = 1;
957              mm[0] = ((double) nlev)/REAL(VECTOR_ELT(Omega, i))[0];              mm[0] = ((double) nlev)/REAL(VECTOR_ELT(Omega, i))[0];
958              mm[1] = F77_CALL(ddot)(&nlev, bi, &ione, bi, &ione);              mm[1] = F77_CALL(ddot)(&nlev, bi, &ione, bi, &ione);
959              mm[2] = F77_CALL(ddot)(&nlev, Di, &ione, Di, &ione);              mm[2] = F77_CALL(ddot)(&nlev, Di, &ione, Di, &ione);
960                mm[3] = 0.
961              for (j = 0; j < p; j++) {              for (j = 0; j < p; j++) {
962                  mm[3] += F77_CALL(ddot)(&RZXrows, RZXi + j * dims[0], &ione,                  mm[3] += F77_CALL(ddot)(&RZXrows, RZXi + j * dims[0], &ione,
963                                          RZXi + j * dims[0], &ione);                                          RZXi + j * dims[0], &ione);
964              }              }
965          } else {          } else {
             double *tmp = Memcpy(Calloc(ncisqr, double),  
                                  REAL(VECTOR_ELT(Omega, i)), ncisqr),  
                 one = 1., zero = 0.;  
   
966              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
967              if (j)              if (j)
968                  error("Omega[[%d]] is not positive definite", i + 1);                  error("Omega[[%d]] is not positive definite", i + 1);
# Line 969  Line 971 
971                  error("Omega[[%d]] is not positive definite", i + 1);                  error("Omega[[%d]] is not positive definite", i + 1);
972              F77_CALL(dsyrk)("U", "N", &nci, &nci, &dlev, tmp, &nci,              F77_CALL(dsyrk)("U", "N", &nci, &nci, &dlev, tmp, &nci,
973                              &zero, mm, &nci);                              &zero, mm, &nci);
             Free(tmp);  
974              mm += ncisqr;       /* \bB_i term */              mm += ncisqr;       /* \bB_i term */
975              F77_CALL(dsyrk)("U", "N", &nci, &nlev, &one, bi, &nci,              F77_CALL(dsyrk)("U", "N", &nci, &nlev, &one, bi, &nci,
976                              &zero, mm, &nci);                              &zero, mm, &nci);
977              mm += ncisqr;     /* Sum of inverses of diagonal blocks */              mm += ncisqr;     /* Sum of inverses of diagonal blocks */
978                /* FIXME: This is wrong.  With more than one grouping
979                 * factor a bVar array will need to be calculated.
980                 * Alternatively, the array could be updated from this. It
981                 * does contain this expression as one factor.*/
982              F77_CALL(dsyrk)("U", "N", &nci, &RZXrows, &one, Di, &nci,              F77_CALL(dsyrk)("U", "N", &nci, &RZXrows, &one, Di, &nci,
983                              &zero, mm, &nci);                              &zero, mm, &nci);
984              mm += ncisqr;       /* Extra term for \vb */              mm += ncisqr;       /* Extra term for \vb */
# Line 984  Line 989 
989              }              }
990          }          }
991      }      }
992        Free(tmp);
993        UNPROTECT(1);
994      return val;      return val;
995  }  }
996    
997  /**  /**
998   * Generate and zero the contents of a list of length nf of arrays of   * Return a length nf list of arrays of dimension (nci, nci, 4).  The
999   * dimension (nci, nci, 4).  The values of these arrays are assigned   * values of these arrays are assigned in lmer_firstDer.
  * in lmer_firstDer.  
1000   *   *
1001   * @param nf number of factors   * @param nf number of factors
1002   * @param nc vector of number of columns per factor   * @param nc vector of number of columns per factor
# Line 1001  Line 1007 
1007  SEXP EM_grad_array(int nf, const int nc[])  SEXP EM_grad_array(int nf, const int nc[])
1008  {  {
1009      SEXP val = PROTECT(allocVector(VECSXP, nf)),      SEXP val = PROTECT(allocVector(VECSXP, nf)),
1010          dd = PROTECT(allocVector(INTSXP, 3));      int i;
     int *dims = INTEGER(dd), i;  
1011    
     dims[2] = 4;  
1012      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
1013          dims[0] = dims[1] = nc[i];          SET_VECTOR_ELT(val, i, alloc3Darray(REALSXP, nc[i], nc[i], 4));
         SET_VECTOR_ELT(val, i, allocArray(REALSXP, dd));  
         memset(REAL(VECTOR_ELT(val, i)), 0, sizeof(double) * nc[i] * nc[i]);  
1014      }      }
1015      UNPROTECT(2);      UNPROTECT(1);
1016      return val;      return val;
1017  }  }
1018    

Legend:
Removed from v.421  
changed lines
  Added in v.422

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