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 565, Tue Feb 22 01:15:45 2005 UTC revision 582, Mon Feb 28 18:15:21 2005 UTC
# Line 157  Line 157 
157          *dims = INTEGER(dd);          *dims = INTEGER(dd);
158    
159      if (length(ncv) != (nf + 1))      if (length(ncv) != (nf + 1))
160          error("length of nc (%d) should be length of facs (%d) + 1",          error(_("length of nc (%d) should be length of facs (%d) + 1"),
161                length(ncv), nf);                length(ncv), nf);
162      SET_VECTOR_ELT(val, 0, NEW_OBJECT(MAKE_CLASS("ssclme")));      SET_VECTOR_ELT(val, 0, NEW_OBJECT(MAKE_CLASS("ssclme")));
163      ssc = VECTOR_ELT(val, 0);      ssc = VECTOR_ELT(val, 0);
# Line 260  Line 260 
260          for (k = 0; k < ncj; k++) {          for (k = 0; k < ncj; k++) {
261              diag = Ap[i + k + 1] - 1;              diag = Ap[i + k + 1] - 1;
262              if (Ai[diag] != i+k)              if (Ai[diag] != i+k)
263                  error("Expected Ai[%d] to be %d (i.e on diagonal) not %d",                  error(_("Expected Ai[%d] to be %d (i.e on diagonal) not %d"),
264                        diag, i+k, Ai[diag]);                        diag, i+k, Ai[diag]);
265              Memcpy(Ax + diag - k, bVj + (i+k-Gpj)*ncj, k + 1);              Memcpy(Ax + diag - k, bVj + (i+k-Gpj)*ncj, k + 1);
266          }          }
# Line 342  Line 342 
342          int *dims = INTEGER(getAttrib(VECTOR_ELT(mmats, i), R_DimSymbol)),          int *dims = INTEGER(getAttrib(VECTOR_ELT(mmats, i), R_DimSymbol)),
343              nci = nc[i];              nci = nc[i];
344          if (nobs != dims[0])          if (nobs != dims[0])
345              error("Expected %d rows in the %d'th model matrix. Got %d",              error(_("Expected %d rows in the %d'th model matrix. Got %d"),
346                    nobs, i+1, dims[0]);                    nobs, i+1, dims[0]);
347          if (nci != dims[1])          if (nci != dims[1])
348              error("Expected %d columns in the %d'th model matrix. Got %d",              error(_("Expected %d columns in the %d'th model matrix. Got %d"),
349                    nci, i+1, dims[1]);                    nci, i+1, dims[1]);
350          Z[i] = REAL(VECTOR_ELT(mmats, i));          Z[i] = REAL(VECTOR_ELT(mmats, i));
351      }      }
# Line 370  Line 370 
370              int fpji = fpj[i] - 1, /* factor indices are 1-based */              int fpji = fpj[i] - 1, /* factor indices are 1-based */
371                  dind = Ap[Gp[j] + fpji * ncj + 1] - 1;                  dind = Ap[Gp[j] + fpji * ncj + 1] - 1;
372              if (Ai[dind] != (Gp[j] + fpji * ncj))              if (Ai[dind] != (Gp[j] + fpji * ncj))
373                  error("logic error in ssclme_update_mm");                  error(_("logic error in ssclme_update_mm"));
374              if (Ncj) {          /* use bVar to accumulate */              if (Ncj) {          /* use bVar to accumulate */
375                  F77_CALL(dsyrk)("U", "T", &ncj, &ione, &one, Zj+i,                  F77_CALL(dsyrk)("U", "T", &ncj, &ione, &one, Zj+i,
376                                  &nobs, &one, bVj + fpji*ncj*ncj, &ncj);                                  &nobs, &one, bVj + fpji*ncj*ncj, &ncj);
# Line 403  Line 403 
403                          break;                          break;
404                      }                      }
405                  }                  }
406                  if (ind < 0) error("logic error in ssclme_update_mm");                  if (ind < 0) error(_("logic error in ssclme_update_mm"));
407                  if (scalar) {   /* update scalars directly */                  if (scalar) {   /* update scalars directly */
408                      Ax[ind] += Zj[i] * Zk[i];                      Ax[ind] += Zj[i] * Zk[i];
409                  } else {                  } else {
# Line 414  Line 414 
414                      for (jj = 0; jj < nck; jj++) {                      for (jj = 0; jj < nck; jj++) {
415                          ind = Apk[fpki * nck + jj] + offset;                          ind = Apk[fpki * nck + jj] + offset;
416                          if (Ai[ind] != row)                          if (Ai[ind] != row)
417                              error("logic error in ssclme_update_mm");                              error(_("logic error in ssclme_update_mm"));
418                          for (ii = 0; ii < ncj; ii++) {                          for (ii = 0; ii < ncj; ii++) {
419                              Ax[ind++] += work[jj * ncj + ii];                              Ax[ind++] += work[jj * ncj + ii];
420                          }                          }
# Line 468  Line 468 
468              for (k = 0; k < ncj; k++) {              for (k = 0; k < ncj; k++) {
469                  diag = Ap[i + k + 1] - 1;                  diag = Ap[i + k + 1] - 1;
470                  if (Ai[diag] != i+k)                  if (Ai[diag] != i+k)
471                      error("Expected Ai[%d] to be %d (i.e on diagonal) not %d",                      error(_("Expected Ai[%d] to be %d (i.e on diagonal) not %d"),
472                            diag, i+k, Ai[diag]);                            diag, i+k, Ai[diag]);
473                  for (ii = 0; ii <= k; ii++) {                  for (ii = 0; ii <= k; ii++) {
474                      xcp[diag + ii - k] += omgj[k*ncj + ii];                      xcp[diag + ii - k] += omgj[k*ncj + ii];
# Line 483  Line 483 
483                      REAL(GET_SLOT(x, Matrix_LxSym)),                      REAL(GET_SLOT(x, Matrix_LxSym)),
484                      D, (int *) NULL, (int *) NULL); /* P & Pinv */                      D, (int *) NULL, (int *) NULL); /* P & Pinv */
485      if (j != n)      if (j != n)
486          error("rank deficiency of ZtZ+W detected at column %d",          error(_("rank deficiency of ZtZ+W detected at column %d"),
487                j + 1);                j + 1);
488      for (j = 0; j < n; j++) DIsqrt[j] = 1./sqrt(D[j]);      for (j = 0; j < n; j++) DIsqrt[j] = 1./sqrt(D[j]);
489      Free(Lnz); Free(Flag); Free(Pattern); Free(Y); Free(xcp);      Free(Lnz); Free(Flag); Free(Pattern); Free(Y); Free(xcp);
# Line 552  Line 552 
552                                          nci * nci),                                          nci * nci),
553                                   &nci, &j);                                   &nci, &j);
554                  if (j)                  if (j)
555                      error("Omega[%d] is not positive definite", i + 1);                      error(_("Omega[%d] is not positive definite"), i + 1);
556                  for (j = 0; j < nci; j++) {                  for (j = 0; j < nci; j++) {
557                      accum += 2 * log(tmp[j * (nci + 1)]);                      accum += 2 * log(tmp[j * (nci + 1)]);
558                  }                  }
# Line 608  Line 608 
608  int ldl_update_ind(int probe, int start, const int ind[])  int ldl_update_ind(int probe, int start, const int ind[])
609  {  {
610      while (ind[start] < probe) start++;      while (ind[start] < probe) start++;
611      if (ind[start] > probe) error("logic error in ldl_inverse");      if (ind[start] > probe) error(_("logic error in ldl_inverse"));
612      return start;      return start;
613  }  }
614    
# Line 696  Line 696 
696                  F77_CALL(dpotrf)("L", &nci, bVi + (j - G1)*nci,                  F77_CALL(dpotrf)("L", &nci, bVi + (j - G1)*nci,
697                                   &nci, &jj);                                   &nci, &jj);
698                  if (jj)         /* should never happen */                  if (jj)         /* should never happen */
699                      error(                      error(_(
700                          "Rank deficient variance matrix at group %d, level %d, error code %d",                          "Rank deficient variance matrix at group %d, level %d, error code %d"),
701                          i + 1, j + 1, jj);                          i + 1, j + 1, jj);
702              }              }
703              Free(tmp); Free(ind);              Free(tmp); Free(ind);
# Line 718  Line 718 
718      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));      int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
719      if (!status[0]) ssclme_factor(x);      if (!status[0]) ssclme_factor(x);
720      if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0]))      if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0]))
721          error("Unable to invert singular factor of downdated X'X");          error(_("Unable to invert singular factor of downdated X'X"));
722      if (!status[1]) {      if (!status[1]) {
723          SEXP          SEXP
724              RZXsl = GET_SLOT(x, Matrix_RZXSym);              RZXsl = GET_SLOT(x, Matrix_RZXSym);
# Line 738  Line 738 
738    
739          F77_CALL(dtrtri)("U", "N", &pp1, RXX, &pp1, &i);          F77_CALL(dtrtri)("U", "N", &pp1, RXX, &pp1, &i);
740          if (i)          if (i)
741              error("DTRTRI returned error code %d", i);              error(_("DTRTRI returned error code %d"), i);
742          F77_CALL(dtrmm)("R", "U", "N", "N", &n, &pp1, &one,          F77_CALL(dtrmm)("R", "U", "N", "N", &n, &pp1, &one,
743                          RXX, &pp1, RZX, &n);                          RXX, &pp1, RZX, &n);
744          for (i = 0; i < pp1; i++) {          for (i = 0; i < pp1; i++) {
# Line 787  Line 787 
787          for (j = Gpi; j < p2; j += nci) {          for (j = Gpi; j < p2; j += nci) {
788              for (k = 0; k < nci; k++) {              for (k = 0; k < nci; k++) {
789                  int jk = j+k, jj = Ap[jk+1] - 1;                  int jk = j+k, jj = Ap[jk+1] - 1;
790                  if (Ai[jj] != jk) error("malformed ZtZ structure");                  if (Ai[jj] != jk) error(_("malformed ZtZ structure"));
791                  mm[k * ncip1] += Ax[jj] * mi;                  mm[k * ncip1] += Ax[jj] * mi;
792              }              }
793          }          }
# Line 936  Line 936 
936                                       REAL(VECTOR_ELT(Omega, i)), ncisq);                                       REAL(VECTOR_ELT(Omega, i)), ncisq);
937                  F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);                  F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
938                  if (j)          /* should never happen */                  if (j)          /* should never happen */
939                      error("DPOTRF returned error code %d on Omega[[%d]]",                      error(_("DPOTRF returned error code %d on Omega[[%d]]"),
940                            j, i+1);                            j, i+1);
941                  for (j = 0; j < nci; j++) {                  for (j = 0; j < nci; j++) {
942                      double diagj = tmp[j * ncip1];                      double diagj = tmp[j * ncip1];
# Line 1001  Line 1001 
1001                                   REAL(VECTOR_ELT(Omega, i)), ncisq);                                   REAL(VECTOR_ELT(Omega, i)), ncisq);
1002              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);              F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
1003              if (j)              /* should never happen */              if (j)              /* should never happen */
1004                  error("DPOTRF returned error code %d on Omega[[%d]]",                  error(_("DPOTRF returned error code %d on Omega[[%d]]"),
1005                        j, i+1);                        j, i+1);
1006              for (j = 0; j < nci; j++) {              for (j = 0; j < nci; j++) {
1007                  double diagj = tmp[j * ncip1];                  double diagj = tmp[j * ncip1];
# Line 1039  Line 1039 
1039      double *cc = REAL(coef);      double *cc = REAL(coef);
1040    
1041      if (length(coef) != coef_length(nf, nc) || !isReal(coef))      if (length(coef) != coef_length(nf, nc) || !isReal(coef))
1042          error("coef must be a numeric vector of length %d",          error(_("coef must be a numeric vector of length %d"),
1043                coef_length(nf, nc));                coef_length(nf, nc));
1044      cind = 0;      cind = 0;
1045      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
# Line 1094  Line 1094 
1094      double *cc = REAL(coef);      double *cc = REAL(coef);
1095    
1096      if (length(coef) != coef_length(nf, nc) || !isReal(coef))      if (length(coef) != coef_length(nf, nc) || !isReal(coef))
1097          error("coef must be a numeric vector of length %d",          error(_("coef must be a numeric vector of length %d"),
1098                coef_length(nf, nc));                coef_length(nf, nc));
1099      cind = 0;      cind = 0;
1100      for (i = 0; i < nf; i++) {      for (i = 0; i < nf; i++) {
# Line 1245  Line 1245 
1245    
1246              F77_CALL(dpotrf)("U", nc + i, vali, nc + i, &info);              F77_CALL(dpotrf)("U", nc + i, vali, nc + i, &info);
1247              if (info)              if (info)
1248                  error("DPOTRF returned error code %d in Omega[[%d]] update",                  error(_("DPOTRF returned error code %d in Omega[[%d]] update"),
1249                        info, i + 1);                        info, i + 1);
1250              F77_CALL(dpotri)("U", nc + i, vali, nc + i, &info);              F77_CALL(dpotri)("U", nc + i, vali, nc + i, &info);
1251              if (info)              if (info)
1252                  error("DPOTRI returned error code %d in Omega[[%d]] update",                  error(_("DPOTRI returned error code %d in Omega[[%d]] update"),
1253                        info, i + 1);                        info, i + 1);
1254          }          }
1255          if (verbose) EMsteps_verbose_print(x, iter + 1, REML);          if (verbose) EMsteps_verbose_print(x, iter + 1, REML);
# Line 1284  Line 1284 
1284    
1285          F77_CALL(dpotrf)("U", &nci, work, &nci, &info);          F77_CALL(dpotrf)("U", &nci, work, &nci, &info);
1286          if (info)          if (info)
1287              error("DPOTRF gave error code %d on Omega[[%d]]", info, i + 1);              error(_("DPOTRF gave error code %d on Omega[[%d]]"), info, i + 1);
1288          F77_CALL(dtrtri)("U", "N", &nci, work, &nci, &info);          F77_CALL(dtrtri)("U", "N", &nci, work, &nci, &info);
1289          if (info)          if (info)
1290              error("DPOTRI gave error code %d on Omega[[%d]]", info, i + 1);              error(_("DPOTRI gave error code %d on Omega[[%d]]"), info, i + 1);
1291          F77_CALL(dsyrk)("U", "N", &nci, &nci, &alpha, work, &nci,          F77_CALL(dsyrk)("U", "N", &nci, &nci, &alpha, work, &nci,
1292                          &beta, REAL(VECTOR_ELT(val, i)), &nci);                          &beta, REAL(VECTOR_ELT(val, i)), &nci);
1293          Free(work);          Free(work);
# Line 1317  Line 1317 
1317    
1318          F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);          F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1319          if (j)          if (j)
1320              error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);              error(_("DPOTRF gave error code %d on Omega[[%d]]"), j, i + 1);
1321          /* calculate and store grad[i,,] %*% t(R) */          /* calculate and store grad[i,,] %*% t(R) */
1322          for (j = 0; j < nci; j++) {          for (j = 0; j < nci; j++) {
1323              F77_CALL(dsymv)("U", &nci, &one, ai, &nci, chol + j, &nci,              F77_CALL(dsymv)("U", &nci, &one, ai, &nci, chol + j, &nci,
# Line 1437  Line 1437 
1437    
1438          F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);          F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
1439          if (j)          if (j)
1440              error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);              error(_("DPOTRF gave error code %d on Omega[[%d]]"), j, i + 1);
1441          Memcpy(tmp, chol, ncisq);          Memcpy(tmp, chol, ncisq);
1442          F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);          F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);
1443          if (j)          if (j)
1444              error("DPOTRI gave error code %d on Omega[[%d]]", j, i + 1);              error(_("DPOTRI gave error code %d on Omega[[%d]]"), j, i + 1);
1445          alpha = (double) -mi;          alpha = (double) -mi;
1446          F77_CALL(dsyrk)("U", "N", &nci, &ki,          F77_CALL(dsyrk)("U", "N", &nci, &ki,
1447                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,                          &one, REAL(VECTOR_ELT(bVar, i)), &nci,
# Line 1565  Line 1565 
1565      double *vv, one = 1.0, zero = 0.0;      double *vv, one = 1.0, zero = 0.0;
1566    
1567      if (nf < 1)      if (nf < 1)
1568          error("null factor list passed to ssclme_fitted");          error(_("null factor list passed to ssclme_fitted"));
1569      nobs = length(VECTOR_ELT(facs, 0));      nobs = length(VECTOR_ELT(facs, 0));
1570      val = PROTECT(allocVector(REALSXP, nobs));      val = PROTECT(allocVector(REALSXP, nobs));
1571      vv = REAL(val);      vv = REAL(val);
# Line 1621  Line 1621 
1621    
1622          F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);          F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);
1623          if (j)                  /* shouldn't happen */          if (j)                  /* shouldn't happen */
1624              error("DPOTRF returned error code %d on Omega[%d]",              error(_("DPOTRF returned error code %d on Omega[%d]"),
1625                    j, i + 1);                    j, i + 1);
1626          F77_CALL(dpotri)("U", &nci, mm, &nci, &j);          F77_CALL(dpotri)("U", &nci, mm, &nci, &j);
1627          if (j)                  /* shouldn't happen */          if (j)                  /* shouldn't happen */
1628              error("DTRTRI returned error code %d on Omega[%d]",              error(_("DTRTRI returned error code %d on Omega[%d]"),
1629                    j, i + 1);                    j, i + 1);
1630          nlme_symmetrize(mm, nci);          nlme_symmetrize(mm, nci);
1631      }      }

Legend:
Removed from v.565  
changed lines
  Added in v.582

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