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 70, Mon Apr 12 12:10:01 2004 UTC revision 100, Sat Apr 17 17:03:59 2004 UTC
# Line 188  Line 188 
188      }      }
189  }  }
190    
 static  
191  void ssclme_fill_LIp(int n, const int Parent[], int LIp[])  void ssclme_fill_LIp(int n, const int Parent[], int LIp[])
192  {  {
193      int *sz = Calloc(n, int), i;      int *sz = Calloc(n, int), i;
# Line 463  Line 462 
462      return R_NilValue;      return R_NilValue;
463  }  }
464    
465  SEXP ssclme_inflate_and_factor(SEXP lme)  SEXP ssclme_inflate_and_factor(SEXP x)
466  {  {
467      SEXP      SEXP
468          GpSlot = GET_SLOT(lme, Matrix_GpSym),          GpSlot = GET_SLOT(x, Matrix_GpSym),
469          Omega = GET_SLOT(lme, Matrix_OmegaSym);          Omega = GET_SLOT(x, Matrix_OmegaSym);
470      int n = INTEGER(GET_SLOT(lme, Matrix_DimSym))[1];      int n = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];
471      int      int
472          *Ai = INTEGER(GET_SLOT(lme, Matrix_iSym)),          *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
473          *Ap = INTEGER(GET_SLOT(lme, Matrix_pSym)),          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
474          *Flag = Calloc(n, int),          *Flag = Calloc(n, int),
475          *Gp = INTEGER(GpSlot),          *Gp = INTEGER(GpSlot),
476          *Lnz = Calloc(n, int),          *Lnz = Calloc(n, int),
477          *Pattern = Calloc(n, int),          *Pattern = Calloc(n, int),
478          *nc = INTEGER(GET_SLOT(lme, Matrix_ncSym)),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
479          j,          j,
480          nf = length(GpSlot) - 1;          nf = length(GpSlot) - 1;
481      double      double
482          *D = REAL(GET_SLOT(lme, Matrix_DSym)),          *D = REAL(GET_SLOT(x, Matrix_DSym)),
483          *DIsqrt = REAL(GET_SLOT(lme, Matrix_DIsqrtSym)),          *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
484          *Y = Calloc(n, double),          *Y = Calloc(n, double),
485          *xcp = Calloc(Ap[n], double);          *xcp = Calloc(Ap[n], double);
486    
487      Memcpy(xcp, REAL(GET_SLOT(lme, Matrix_xSym)), Ap[n]);      Memcpy(xcp, REAL(GET_SLOT(x, Matrix_xSym)), Ap[n]);
488      for (j = 0; j < nf; j++) {      for (j = 0; j < nf; j++) {
489          int  diag, i, ii, k, G2 = Gp[j + 1], ncj = nc[j];          int  diag, i, ii, k, G2 = Gp[j + 1], ncj = nc[j];
490          double *omgj = REAL(VECTOR_ELT(Omega, j));          double *omgj = REAL(VECTOR_ELT(Omega, j));
# Line 503  Line 502 
502          }          }
503      }      }
504      j = ldl_numeric(n, Ap, Ai, xcp,      j = ldl_numeric(n, Ap, Ai, xcp,
505                      INTEGER(GET_SLOT(lme, Matrix_LpSym)),                      INTEGER(GET_SLOT(x, Matrix_LpSym)),
506                      INTEGER(GET_SLOT(lme, Matrix_ParentSym)),                      INTEGER(GET_SLOT(x, Matrix_ParentSym)),
507                      Lnz, INTEGER(GET_SLOT(lme, Matrix_LiSym)),                      Lnz, INTEGER(GET_SLOT(x, Matrix_LiSym)),
508                      REAL(GET_SLOT(lme, Matrix_LxSym)),                      REAL(GET_SLOT(x, Matrix_LxSym)),
509                      D, Y, Pattern, Flag,                      D, Y, Pattern, Flag,
510                      (int *) NULL, (int *) NULL); /* P & Pinv */                      (int *) NULL, (int *) NULL); /* P & Pinv */
511      if (j != n)      if (j != n)
# Line 517  Line 516 
516      return R_NilValue;      return R_NilValue;
517  }  }
518    
519  SEXP ssclme_factor(SEXP lme)  SEXP ssclme_factor(SEXP x)
520  {  {
521      int *status = LOGICAL(GET_SLOT(lme, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
522    
523      if (!status[0]) {      if (!status[0]) {
524          SEXP          SEXP
525              GpSlot = GET_SLOT(lme, Matrix_GpSym),              GpSlot = GET_SLOT(x, Matrix_GpSym),
526              Omega = GET_SLOT(lme, Matrix_OmegaSym);              Omega = GET_SLOT(x, Matrix_OmegaSym);
527          int          int
528              *Gp = INTEGER(GpSlot),              *Gp = INTEGER(GpSlot),
529              *Li = INTEGER(GET_SLOT(lme, Matrix_LiSym)),              *Li = INTEGER(GET_SLOT(x, Matrix_LiSym)),
530              *Lp = INTEGER(GET_SLOT(lme, Matrix_LpSym)),              *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
531              *nc = INTEGER(GET_SLOT(lme, Matrix_ncSym)),              *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
532              i,              i,
533              n = INTEGER(GET_SLOT(lme, Matrix_DimSym))[1],              n = INTEGER(GET_SLOT(x, Matrix_DimSym))[1],
534              nf = length(GpSlot) - 1,              nf = length(GpSlot) - 1,
535              nobs = nc[nf + 1],              nobs = nc[nf + 1],
536              nreml = nobs + 1 - nc[nf],              nreml = nobs + 1 - nc[nf],
537              pp1 = nc[nf],              pp1 = nc[nf],
538              pp2 = pp1 + 1;              pp2 = pp1 + 1;
539          double          double
540              *D = REAL(GET_SLOT(lme, Matrix_DSym)),              *D = REAL(GET_SLOT(x, Matrix_DSym)),
541              *DIsqrt = REAL(GET_SLOT(lme, Matrix_DIsqrtSym)),              *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
542              *Lx = REAL(GET_SLOT(lme, Matrix_LxSym)),              *Lx = REAL(GET_SLOT(x, Matrix_LxSym)),
543              *RXX = REAL(GET_SLOT(lme, Matrix_RXXSym)),              *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
544              *RZX = REAL(GET_SLOT(lme, Matrix_RZXSym)),              *RZX = REAL(GET_SLOT(x, Matrix_RZXSym)),
545              *dcmp = REAL(getAttrib(lme, Matrix_devCompSym)),              *dcmp = REAL(getAttrib(x, Matrix_devCompSym)),
546              *deviance = REAL(getAttrib(lme, Matrix_devianceSym)),              *deviance = REAL(getAttrib(x, Matrix_devianceSym)),
547              minus1 = -1.,              minus1 = -1.,
548              one = 1.;              one = 1.;
549    
550          ssclme_inflate_and_factor(lme);          ssclme_inflate_and_factor(x);
551                                  /* Accumulate logdet of ZtZ+W */                                  /* Accumulate logdet of ZtZ+W */
552          dcmp[0] = dcmp[1] = dcmp[2] = dcmp[3] = 0.;          dcmp[0] = dcmp[1] = dcmp[2] = dcmp[3] = 0.;
553          for (i = 0; i < n; i++) dcmp[0] += log(D[i]);          for (i = 0; i < n; i++) dcmp[0] += log(D[i]);
# Line 578  Line 577 
577              }              }
578          }          }
579                                  /* ldl_lsolve on Z'X */                                  /* ldl_lsolve on Z'X */
580          Memcpy(RZX, REAL(GET_SLOT(lme, Matrix_ZtXSym)), n * pp1);          Memcpy(RZX, REAL(GET_SLOT(x, Matrix_ZtXSym)), n * pp1);
581          for (i = 0; i < pp1; i++) {          for (i = 0; i < pp1; i++) {
582              int j;              int j;
583              double *RZXi = RZX + i * n;              double *RZXi = RZX + i * n;
# Line 586  Line 585 
585              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
586          }          }
587                                  /* downdate and factor X'X */                                  /* downdate and factor X'X */
588          Memcpy(RXX, REAL(GET_SLOT(lme, Matrix_XtXSym)), pp1 * pp1);          Memcpy(RXX, REAL(GET_SLOT(x, Matrix_XtXSym)), pp1 * pp1);
589          F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,          F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,
590                          RZX, &n, &one, RXX, &pp1);                          RZX, &n, &one, RXX, &pp1);
591          F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);          F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);
# Line 624  Line 623 
623   * @return R_NilValue (x is updated in place)   * @return R_NilValue (x is updated in place)
624    
625   */   */
626    static
627  SEXP ldl_inverse(SEXP x)  SEXP ldl_inverse(SEXP x)
628  {  {
629      SEXP      SEXP
# Line 791  Line 791 
791      }      }
792      return R_NilValue;      return R_NilValue;
793  }  }
794    SEXP ssclme_invert(SEXP x)
 SEXP ssclme_invert(SEXP lme)  
795  {  {
796      int *status = LOGICAL(GET_SLOT(lme, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
797      if (!status[0]) ssclme_factor(lme);      if (!status[0]) ssclme_factor(x);
798      if (!status[1]) {      if (!status[1]) {
799          SEXP          SEXP
800              RZXsl = GET_SLOT(lme, Matrix_RZXSym);              RZXsl = GET_SLOT(x, Matrix_RZXSym);
801          int          int
802              *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),              *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
803              *Li = INTEGER(GET_SLOT(lme, Matrix_LiSym)),              *Li = INTEGER(GET_SLOT(x, Matrix_LiSym)),
804              *Lp = INTEGER(GET_SLOT(lme, Matrix_LpSym)),              *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
805              i,              i,
806              n = dims[0],              n = dims[0],
807              pp1 = dims[1];              pp1 = dims[1];
808          double          double
809              *DIsqrt = REAL(GET_SLOT(lme, Matrix_DIsqrtSym)),              *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
810              *Lx = REAL(GET_SLOT(lme, Matrix_LxSym)),              *Lx = REAL(GET_SLOT(x, Matrix_LxSym)),
811              *RXX = REAL(GET_SLOT(lme, Matrix_RXXSym)),              *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
812              *RZX = REAL(RZXsl),              *RZX = REAL(RZXsl),
813              one = 1.;              one = 1.;
814    
# Line 823  Line 822 
822              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];              for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
823              ldl_ltsolve(n, RZXi, Lp, Li, Lx);              ldl_ltsolve(n, RZXi, Lp, Li, Lx);
824          }          }
825          ldl_inverse(lme);          ldl_inverse(x);
826          status[1] = 1;          status[1] = 1;
827      }      }
828      return R_NilValue;      return R_NilValue;
# Line 1192  Line 1191 
1191                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,                              &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,
1192                              &one, vali, &nci);                              &one, vali, &nci);
1193              if (REML) {              if (REML) {
1194                  int mp = mi * p;                  int j;
1195                  F77_CALL(dsyrk)("U", "N", &nci, &mp,                  for (j = 0; j < p; j++) {
1196                                  &alpha, RZX + Gp[i], &nci,                      F77_CALL(dsyrk)("U", "N", &nci, &mi,
1197                                    &alpha, RZX + Gp[i] + j*n, &nci,
1198                                  &one, vali, &nci);                                  &one, vali, &nci);
1199              }              }
1200                }
1201              F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);              F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);
1202              if (info) error("DPOTRF returned error code %d", info);              if (info) error("DPOTRF returned error code %d", info);
1203              F77_CALL(dpotri)("U", &nci, vali, &nci, &info);              F77_CALL(dpotri)("U", &nci, vali, &nci, &info);
1204              if (info) error("DPOTRF returned error code %d", info);              if (info) error("DPOTRI returned error code %d", info);
1205          }          }
1206          status[0] = status[1] = 0;          status[0] = status[1] = 0;
1207      }      }
# Line 1208  Line 1209 
1209      return R_NilValue;      return R_NilValue;
1210  }  }
1211    
1212    SEXP ssclme_gradient(SEXP x, SEXP REMLp)
1213    {
1214        SEXP
1215            RZXsl = GET_SLOT(x, Matrix_RZXSym),
1216            ans = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym))),
1217            ncsl = GET_SLOT(x, Matrix_ncSym),
1218            bVar = GET_SLOT(x, Matrix_bVarSym);
1219        int
1220            *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1221            *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1222            *nc = INTEGER(ncsl),
1223            REML = asLogical(REMLp),
1224            i, info,
1225            n = dims[0],
1226            nf = length(ncsl) - 2,
1227            nobs = nc[nf + 1],
1228            p,
1229            pp1 = dims[1];
1230        double
1231            *RZX = REAL(RZXsl),
1232            *b,
1233            alpha,
1234            one = 1.;
1235    
1236        p = pp1 - 1;
1237        b = RZX + p * n;
1238        ssclme_invert(x);
1239        for (i = 0; i < nf; i++) {
1240            int ki = Gp[i+1] - Gp[i],
1241                nci = nc[i],
1242                mi = ki/nci;
1243            double
1244                *vali = REAL(VECTOR_ELT(ans, i));
1245    
1246            F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);
1247            if (info)
1248                error("DPOTRF returned error code %d for component %d of Omega",
1249                      info, i + 1);
1250            F77_CALL(dpotri)("U", &nci, vali, &nci, &info);
1251            if (info)
1252                error("DPOTRI returned error code %d for component %d of Omega",
1253                      info, i + 1);
1254            alpha = (double) -mi;
1255            F77_CALL(dsyrk)("U", "N", &nci, &ki,
1256                            &one, REAL(VECTOR_ELT(bVar, i)), &nci,
1257                            &alpha, vali, &nci);
1258            alpha = ((double)(REML?(nobs-p):nobs));
1259            F77_CALL(dsyrk)("U", "N", &nci, &mi,
1260                            &alpha, b + Gp[i], &nci,
1261                            &one, vali, &nci);
1262            if (REML) {
1263                int j;
1264                for (j = 0; j < p; j++) {
1265                    F77_CALL(dsyrk)("U", "N", &nci, &mi,
1266                                    &one, RZX + Gp[i] + j*n, &nci,
1267                                    &one, vali, &nci);
1268                }
1269            }
1270        }
1271        UNPROTECT(1);
1272        return ans;
1273    }
1274    
1275  SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats)  SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats)
1276  {  {
1277      SEXP val, b;      SEXP val, b;
# Line 1286  Line 1350 
1350      return val;      return val;
1351  }  }
1352    
 SEXP ssclme_L_LI_sizes(SEXP facs)  
 {  
     SEXP ans = PROTECT(allocVector(INTSXP, 4)),  
         ctab;  
     int *Ai, *Ap, *Lp, *Parent, *aa = INTEGER(ans), *dims, *perm,  
         nzcol;  
                                 /* Pairwise cross-tabulation */  
     ctab = PROTECT(sscCrosstab(facs, ScalarLogical(1)));  
     Ai = INTEGER(GET_SLOT(ctab, Matrix_iSym));  
     Ap = INTEGER(GET_SLOT(ctab, Matrix_pSym));  
     dims = INTEGER(GET_SLOT(ctab, Matrix_DimSym));  
     nzcol = dims[1];  
     Lp = Calloc(nzcol + 1, int);  
     Parent = Calloc(nzcol, int);  
     ldl_symbolic(nzcol, Ap, Ai, Lp, Parent,  
                  (int *) R_alloc(nzcol, sizeof(int)), /* Lnz */  
                  (int *) R_alloc(nzcol, sizeof(int)), /* Flag */  
                  (int *) NULL, (int *) NULL); /* P & Pinv */  
     aa[0] = Lp[nzcol];  
     ssclme_fill_LIp(nzcol, Parent, Lp);  
     aa[1] = Lp[nzcol];  
     perm = Calloc(nzcol, int);  
     ssc_metis_order(nzcol, Ap, Ai, perm, Parent);  
     ssc_symbolic_permute(nzcol, 1, Parent, Ap, Ai);  
     ldl_symbolic(nzcol, Ap, Ai, Lp, Parent,  
                  (int *) R_alloc(nzcol, sizeof(int)), /* Lnz */  
                  (int *) R_alloc(nzcol, sizeof(int)), /* Flag */  
                  (int *) NULL, (int *) NULL); /* P & Pinv */  
     aa[2] = Lp[nzcol];  
     ssclme_fill_LIp(nzcol, Parent, Lp);  
     aa[3] = Lp[nzcol];  
     Free(perm); Free(Parent); Free(Lp);  
     UNPROTECT(2);  
     return ans;  
 }  

Legend:
Removed from v.70  
changed lines
  Added in v.100

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