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 76, Wed Apr 14 20:43:40 2004 UTC revision 97, Sat Apr 17 16:57:42 2004 UTC
# Line 462  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 502  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 516  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 577  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 585  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 623  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 790  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 822  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;

Legend:
Removed from v.76  
changed lines
  Added in v.97

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