SCM

SCM Repository

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

Diff of /pkg/src/lmeRep.c

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

revision 372, Sun Dec 5 13:53:29 2004 UTC revision 373, Sun Dec 5 13:54:32 2004 UTC
# Line 1  Line 1 
1  #include "lmeRep.h"  #include "lmeRep.h"
2    
3  /**  /**
4   * Calculate the length of the parameter vector, which is called coef   * Calculate the length of the parameter vector (historically called "coef"
5   * for historical reasons.   * even though these are not coefficients).
6   *   *
7   * @param nf number of factors   * @param nf number of factors
8   * @param nc number of columns in the model matrices for each factor   * @param nc number of columns in the model matrices for each factor
# Line 624  Line 624 
624  }  }
625    
626  /**  /**
627   * Copy the diagonal blocks from ZZx to D and inflate with Omega.  Update   * Copy the diagonal blocks from ZZx to D, inflate by Omega and factor.
628   * dcmp[1] in the process.   * Update dcmp[1] in the process.
629   *   *
630   * @param nf number of grouping factors   * @param nf number of grouping factors
631   * @param ZZxP pointer to the ZZx list   * @param ZZxP pointer to the ZZx list
# Line 650  Line 650 
650          if (nci == 1) {          if (nci == 1) {
651              dcmp[1] += nlev * log(Omega[0]);              dcmp[1] += nlev * log(Omega[0]);
652              for (j = 0; j < nlev; j++) { /* inflate and factor */              for (j = 0; j < nlev; j++) { /* inflate and factor */
653                  D[j] = sqrt(ZZx[j] + Omega[0]);                  D[j] = ZZx[j] + Omega[0];
654              }              }
655          } else {          } else {
656              double *tmp = Memcpy(Calloc(ncisqr, double), Omega, ncisqr);              double *tmp = Memcpy(Calloc(ncisqr, double), Omega, ncisqr);
# Line 671  Line 671 
671                      *ZZxj = ZZx + j * ncisqr;                      *ZZxj = ZZx + j * ncisqr;
672    
673                  for (jj = 0; jj < nci; jj++) { /* Copy ZZx to Dj and inflate */                  for (jj = 0; jj < nci; jj++) { /* Copy ZZx to Dj and inflate */
674                      for (ii = 0; ii <= jj; ii++) {                      for (ii = 0; ii < nci; ii++) { /* upper triangle only */
675                          Dj[ii + jj * nci]                          Dj[ii + jj * nci] = (ii > jj) ? 0. :
676                              = ZZxj[ii + jj * nci] + Omega[ii + jj * nci];                              ZZxj[ii + jj * nci] + Omega[ii + jj * nci];
677                      }                      }
                     for (ii = jj + 1; ii < nci; ii++)  
                         Dj[ii + jj * nci] = 0.;  
678                  }                  }
679              }              }
680          }          }
# Line 839  Line 837 
837    
838              if (i) update_D_L(i, nc, LP, DP);              if (i) update_D_L(i, nc, LP, DP);
839              if (nci == 1) {              if (nci == 1) {
840                  for (j = 0; j < nlev; j++)                  for (j = 0; j < nlev; j++) {
841                        D[j] = sqrt(D[j]);
842                      dcmp[0] += 2. * log(D[j]);                      dcmp[0] += 2. * log(D[j]);
843                    }
844              } else {              } else {
845                  for (j = 0; j < nlev; j++) {                  for (j = 0; j < nlev; j++) {
846                      double *Dj = D + j * ncisqr;                      double *Dj = D + j * ncisqr;
# Line 878  Line 878 
878          F77_CALL(dpotrf)("U", dims + 1, RXX, dims + 1, &j);          F77_CALL(dpotrf)("U", dims + 1, RXX, dims + 1, &j);
879          if (j) {          if (j) {
880              warning("Leading minor of size %d of downdated X'X is indefinite",              warning("Leading minor of size %d of downdated X'X is indefinite",
881                      j + 1);                      j);
882              dcmp[2] = dcmp[3] = deviance[0] = deviance[1] = NA_REAL;              dcmp[2] = dcmp[3] = deviance[0] = deviance[1] = NA_REAL;
883          } else {          } else {
884              for (j = 0; j < (dims[1] - 1); j++) /* 2 logDet(RXX) */              for (j = 0; j < (dims[1] - 1); j++) /* 2 logDet(RXX) */

Legend:
Removed from v.372  
changed lines
  Added in v.373

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