# SCM Repository

[matrix] Diff of /pkg/src/ssclme.c
 [matrix] / pkg / src / ssclme.c

# Diff of /pkg/src/ssclme.c

revision 158, Sun May 9 17:06:12 2004 UTC revision 159, Sun May 9 22:06:53 2004 UTC
# Line 99  Line 99
99      }      }
100  }  }
101
102  void ssclme_fill_LIp(int n, const int Parent[], int LIp[])  /**
103     * Calculate and store the maximum number of off-diagonal elements in
104     * the inverse of L, based on the elimination tree.  The maximum is
105     * itself stored in the Parent array.  (FIXME: come up with a better design.)
106     *
107     * @param n number of columns in the matrix
108     * @param Parent elimination tree for the matrix
109     */
110    static void ssclme_calc_maxod(int n, int Parent[])
111  {  {
112      int *sz = Calloc(n, int), i;      int *sz = Calloc(n, int), i, mm = -1;
113      for (i = n - 1; i >= 0; i--) {      for (i = n - 1; i >= 0; i--) {
114          sz[i] = (Parent[i] < 0) ? 0 : 1 + sz[Parent[i]];          sz[i] = (Parent[i] < 0) ? 0 : (1 + sz[Parent[i]]);
115            if (sz[i] > mm) mm = sz[i];
116      }      }
117      LIp[0] = 0;      Parent[n] = mm;
for (i = 0; i < n; i++) LIp[i+1] = LIp[i] + sz[i];
118      Free(sz);      Free(sz);
119  }  }
120
static
void ssclme_fill_LIi(int n, const int Parent[], const int LIp[], int LIi[])
{
int i;
for (i = n; i > 0; i--) {
int im1 = i - 1, Par = Parent[im1];
if (Par >= 0) {
LIi[LIp[im1]] = Par;
Memcpy(LIi + LIp[im1] + 1, LIi + LIp[Par],
LIp[Par + 1] - LIp[Par]);
}
}
}

121  SEXP  SEXP
122  ssclme_create(SEXP facs, SEXP ncv, SEXP threshold)  ssclme_create(SEXP facs, SEXP ncv, SEXP threshold)
123  {  {
124      SEXP ctab, nms, ssc, tmp,      SEXP ctab, nms, ssc, tmp,
125          val = PROTECT(allocVector(VECSXP, 2)),          val = PROTECT(allocVector(VECSXP, 2)),
126          dd = PROTECT(allocVector(INTSXP, 3));   /* dimensions of 3-D arrays */          dd = PROTECT(allocVector(INTSXP, 3));   /* dimensions of 3-D arrays */
127      int *Ai, *Ap, *Gp, *LIp, *Lp, *Parent,      int *Ai, *Ap, *Gp, *Lp, *Parent,
128          *nc, Lnz, i, nf = length(facs), nzcol, pp1,          *nc, Lnz, i, nf = length(facs), nzcol, pp1,
129          *dims = INTEGER(dd);          *dims = INTEGER(dd);
130
# Line 170  Line 164
164             sizeof(double) * pp1 * pp1);             sizeof(double) * pp1 * pp1);
165      SET_SLOT(ssc, Matrix_LpSym, allocVector(INTSXP, nzcol + 1));      SET_SLOT(ssc, Matrix_LpSym, allocVector(INTSXP, nzcol + 1));
166      Lp = INTEGER(GET_SLOT(ssc, Matrix_LpSym));      Lp = INTEGER(GET_SLOT(ssc, Matrix_LpSym));
167      SET_SLOT(ssc, Matrix_ParentSym, allocVector(INTSXP, nzcol));      SET_SLOT(ssc, Matrix_ParentSym, allocVector(INTSXP, nzcol + 1));
168      Parent = INTEGER(GET_SLOT(ssc, Matrix_ParentSym));      Parent = INTEGER(GET_SLOT(ssc, Matrix_ParentSym));
169      SET_SLOT(ssc, Matrix_DSym, allocVector(REALSXP, nzcol));      SET_SLOT(ssc, Matrix_DSym, allocVector(REALSXP, nzcol));
170      SET_SLOT(ssc, Matrix_DIsqrtSym, allocVector(REALSXP, nzcol));      SET_SLOT(ssc, Matrix_DIsqrtSym, allocVector(REALSXP, nzcol));
# Line 178  Line 172
172                   (int *) R_alloc(nzcol, sizeof(int)), /* Lnz */                   (int *) R_alloc(nzcol, sizeof(int)), /* Lnz */
173                   (int *) R_alloc(nzcol, sizeof(int)), /* Flag */                   (int *) R_alloc(nzcol, sizeof(int)), /* Flag */
174                   (int *) NULL, (int *) NULL); /* P & Pinv */                   (int *) NULL, (int *) NULL); /* P & Pinv */
175        ssclme_calc_maxod(nzcol, Parent);
176      Lnz = Lp[nzcol];      Lnz = Lp[nzcol];
177      SET_SLOT(ssc, Matrix_LiSym, allocVector(INTSXP, Lnz));      SET_SLOT(ssc, Matrix_LiSym, allocVector(INTSXP, Lnz));
178      SET_SLOT(ssc, Matrix_LxSym, allocVector(REALSXP, Lnz));      SET_SLOT(ssc, Matrix_LxSym, allocVector(REALSXP, Lnz));
# Line 215  Line 210
210          memset(REAL(VECTOR_ELT(tmp, i)), 0,          memset(REAL(VECTOR_ELT(tmp, i)), 0,
211                 sizeof(double) * nci * nci * mi);                 sizeof(double) * nci * nci * mi);
212      }      }
SET_SLOT(ssc, Matrix_LIpSym, allocVector(INTSXP, nzcol + 1));
LIp = INTEGER(GET_SLOT(ssc, Matrix_LIpSym));
ssclme_fill_LIp(nzcol, Parent, LIp);
if (asInteger(threshold) > (Lnz = LIp[nzcol])) {
SET_SLOT(ssc, Matrix_LIiSym, allocVector(INTSXP, Lnz));
ssclme_fill_LIi(nzcol, Parent, LIp,
INTEGER(GET_SLOT(ssc, Matrix_LIiSym)));
SET_SLOT(ssc, Matrix_LIxSym, allocVector(REALSXP, Lnz));
memset(REAL(GET_SLOT(ssc, Matrix_LIxSym)), 0,
sizeof(double) * Lnz);
}
213      UNPROTECT(2);      UNPROTECT(2);
214      return val;      return val;
215  }  }
# Line 543  Line 527
527  }  }
528
529  /**  /**
530   * Create the inverse of L and update the diagonal blocks of the inverse   * Update the diagonal blocks of the inverse of LDL' (=Z'Z+W)
* of LDL' (=Z'Z+W)
531   *   *
532   * @param x pointer to an ssclme object   * @param x pointer to an ssclme object
533   *   *
# Line 556  Line 539
539  {  {
540      SEXP      SEXP
541          Gpsl = GET_SLOT(x, Matrix_GpSym),          Gpsl = GET_SLOT(x, Matrix_GpSym),
LIisl = GET_SLOT(x, Matrix_LIiSym),
LIpsl = GET_SLOT(x, Matrix_LIpSym),
542          bVar = GET_SLOT(x, Matrix_bVarSym);          bVar = GET_SLOT(x, Matrix_bVarSym);
543      int *Gp = INTEGER(Gpsl),      int *Gp = INTEGER(Gpsl),
544          *Li,          *Parent = INTEGER(GET_SLOT(x, Matrix_ParentSym)),
545          *LIp = INTEGER(LIpsl), *Lp,          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
546          i,          i,
547          nf = length(Gpsl) - 1,          nf = length(Gpsl) - 1,
548          nzc = length(LIpsl) - 1;          nzc = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];
549      double      int maxod = Parent[nzc];
550          *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),      double *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym));
*Lx;
551
552      ssclme_factor(x);      ssclme_factor(x);
553      if (LIp[nzc] == 0) {        /* L and LI are the identity */      if (maxod == 0) {           /* L and L^{-1} are the identity */
554          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
555              Memcpy(REAL(VECTOR_ELT(bVar, i)), DIsqrt + Gp[i],              Memcpy(REAL(VECTOR_ELT(bVar, i)), DIsqrt + Gp[i],
556                     Gp[i+1] - Gp[i]);                     Gp[i+1] - Gp[i]);
557          }          }
558          return R_NilValue;      } else {
559      }          int *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
560      Lp = INTEGER(GET_SLOT(x, Matrix_LpSym));              *Li = INTEGER(GET_SLOT(x, Matrix_LiSym));
Li = INTEGER(GET_SLOT(x, Matrix_LiSym));
Lx = REAL(GET_SLOT(x, Matrix_LxSym));
if (length(LIisl) == LIp[nzc]) { /* LIi is filled */
int *LIi = INTEGER(LIisl),
*nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
j, jj, k, kk, p1, p2, pi1, pi2;
561
562          double *LIx = REAL(GET_SLOT(x, Matrix_LIxSym)),          double one = 1.0, zero = 0.,
563              one = 1., zero = 0.;              *Lx = REAL(GET_SLOT(x, Matrix_LxSym));

memset(LIx, 0, sizeof(double) * LIp[nzc]);
/* calculate inverse */
for (i = 0; i < nzc; i++) {
p1 = Lp[i]; p2 = Lp[i+1]; pi1 = LIp[i]; pi2 = LIp[i+1];
/* initialize from unit diagonal term */
kk = pi1;
for (j = p1; j < p2; j++) {
k = Li[j];
while (LIi[kk] < k && kk < pi2) kk++;
if (LIi[kk] != k) error("logic error in ldl_inverse");
LIx[kk] = -Lx[j];
}
for (j = pi1; j < pi2; j++) {
jj = LIi[j];
p1 = Lp[jj]; p2 = Lp[jj+1];
kk = j;
for (jj = p1; jj < p2; jj++) {
k = Li[jj];
while (LIi[kk] < k && kk < pi2) kk++;
if (LIi[kk] != k) error("logic error in ldl_inverse");
LIx[kk] -= Lx[jj]*LIx[j];
}
}
}
for (i = 0; i < nf; i++) { /* accumulate bVar */
int G1 = Gp[i], G2 = Gp[i+1], j, k, kk,
nci = nc[i], nr, nr1, rr;
double *bVi = REAL(VECTOR_ELT(bVar, i)), *tmp;

nr = -1;
for (j = G1; j < G2; j += nci) {
rr = 1 + LIp[j + 1] - LIp[j];
if (rr > nr) nr = rr;
}
tmp = Calloc(nr * nci, double); /* scratch storage */
nr1 = nr + 1;
/* initialize bVi to zero */
memset(bVi, 0, sizeof(double) * (G2 - G1) * nci);
for (j = G1; j < G2; j += nci) {
memset(tmp, 0, sizeof(double) * nr * nci);
rr = 1 + LIp[j + 1] - LIp[j];
for (k = 0; k < nci; k++) { /* copy columns */
tmp[k * nr1] = 1.; /* (unstored) diagonal elt  */
Memcpy(tmp + k*nr1 + 1, LIx + LIp[j + k], rr - k - 1);
}
/* scale the rows */
tmp[0] = DIsqrt[j]; /* first row only has one non-zero */
for (kk = 1; kk < rr; kk++) {
for (k = 0; k < nci; k++) {
tmp[k * nr + kk] *= DIsqrt[LIi[LIp[j] + kk - 1]];
}
}
F77_CALL(dsyrk)("L", "T", &nci, &rr, &one, tmp, &nr,
&zero, bVi + (j - G1) * nci, &nci);
F77_CALL(dpotrf)("L", &nci, bVi + (j - G1) * nci,
&nci, &kk);
if (kk)         /* should never happen */
error(
"Rank deficient variance matrix at group %d, level %d",
i + 1, j + 1);
}
Free(tmp);
}
return R_NilValue;
}
if (length(LIisl)) error("logic error in ssclme_ldl_inverse");
else {                  /* LIi and LIx are too big and not used */
int info, maxod;
int *Parent = INTEGER(GET_SLOT(x, Matrix_ParentSym));
int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym));
double one = 1.0, zero = 0.;

maxod = -1;         /* determine maximum # of off-diagonals */
for (i = 0; i < nzc; i++) { /* in a column of L^{-1} */
int ci = LIp[i+1] - LIp[i];
if (ci > maxod) maxod = ci;
}
/* FIXME: Don't bother storing LIp, store maxod only */
564
565          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
566              int j, jj, k, kk, nci = nc[i], nr, p, p2, pj, pp,              int j, jj, k, kk, nci = nc[i], nr, p, p2, pj, pp,
567                  m = maxod + 1,                  m = maxod + 1,
568                  *ind = Calloc(m, int);                  *ind = Calloc(m, int), G1 = Gp[i], G2 = Gp[i+1];
569              double              double
570                  *tmp = Calloc(m * nci, double),                  *tmp = Calloc(m * nci, double),
571                  *mpt = REAL(VECTOR_ELT(bVar, i));                  *bVi = REAL(VECTOR_ELT(bVar, i));
572
573              for (j = Gp[i]; j < Gp[i+1]; j += nci) {                                  /* initialize bVi to zero */
574                memset(bVi, 0, sizeof(double) * (G2 - G1) * nci);
575
576                for (j = G1; j < G2; j += nci) {
577                  kk = 0;         /* ind gets indices of non-zeros */                  kk = 0;         /* ind gets indices of non-zeros */
578                  jj = j;         /* in this block of columns */                  jj = j;         /* in this block of columns */
579                  while (jj >= 0) {                  while (jj >= 0) {
# Line 706  Line 604
604                      }                      }
605                  }                  }
606                  F77_CALL(dsyrk)("L", "T", &nci, &nr, &one, tmp, &nr,                  F77_CALL(dsyrk)("L", "T", &nci, &nr, &one, tmp, &nr,
607                                  &zero, mpt + (j - Gp[i])*nci, &nci);                                  &zero, bVi + (j - G1)*nci, &nci);
608                  F77_CALL(dpotrf)("L", &nci, mpt + (j - Gp[i])*nci,                  F77_CALL(dpotrf)("L", &nci, bVi + (j - G1)*nci,
609                                   &nci, &info);                                   &nci, &jj);
610                  if (info)       /* should never happen */                  if (jj)         /* should never happen */
611                      error(                      error(
612                          "Rank deficient variance matrix at group %d, level %d",                          "Rank deficient variance matrix at group %d, level %d, error code %d",
613                          i + 1, j + 1);                          i + 1, j + 1, jj);
614              }              }
615              Free(tmp); Free(ind);              Free(tmp); Free(ind);
616          }          }
# Line 1343  Line 1241
1241          Dim = GET_SLOT(x, Matrix_DimSym);          Dim = GET_SLOT(x, Matrix_DimSym);
1242      int i, nf = length(Omega), nz = INTEGER(Dim)[1];      int i, nf = length(Omega), nz = INTEGER(Dim)[1];
1243      SEXP copy[] = {Matrix_DSym, Matrix_DIsqrtSym, Matrix_DimSym,      SEXP copy[] = {Matrix_DSym, Matrix_DIsqrtSym, Matrix_DimSym,
1244                     Matrix_GpSym, Matrix_LIiSym, Matrix_LIpSym,                     Matrix_GpSym, Matrix_LiSym, Matrix_LpSym,
Matrix_LIxSym, Matrix_LiSym, Matrix_LpSym,
1245                     Matrix_LxSym, Matrix_OmegaSym, Matrix_ParentSym,                     Matrix_LxSym, Matrix_OmegaSym, Matrix_ParentSym,
Matrix_LIxSym, Matrix_LiSym, Matrix_LpSym,
1246                     Matrix_bVarSym, Matrix_devianceSym,                     Matrix_bVarSym, Matrix_devianceSym,
1247                     Matrix_devCompSym, Matrix_iSym, Matrix_ncSym,                     Matrix_devCompSym, Matrix_iSym, Matrix_ncSym,
1248                     Matrix_statusSym, Matrix_pSym, Matrix_xSym};                     Matrix_statusSym, Matrix_pSym, Matrix_xSym};
1249
1250      for (i = 0; i < 23; i++)      for (i = 0; i < 17; i++)
1251          SET_SLOT(ans, copy[i], duplicate(GET_SLOT(x, copy[i])));          SET_SLOT(ans, copy[i], duplicate(GET_SLOT(x, copy[i])));
1252
1253      INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;      INTEGER(GET_SLOT(ans, Matrix_ncSym))[nf] = 1;

Legend:
 Removed from v.158 changed lines Added in v.159