#include "ssclme.h"

/** 
 * Check for a nested series of grouping factors in the sparse,
 *  symmetric representation of the pairwise cross-tabulations.
 * 
 * @param n size of pairwise cross-tabulation matrix
 * @param nf number of groups of columns in pairwise cross-tabulation
 * @param upper non-zero if the upper triangle is stored
 * @param Ap array of pointers to columns
 * @param Ai row indices
 * @param Gp array of pointers to groups
 * 
 * @return 0 for non-nested groups, 1 for nested groups
 */
static
int ctab_isNested(int n, int nf, int upper,
		  const int Ap[], const int Ai[], const int Gp[])
{
    if (nf > 1) {  /* single factor always nested */
	int  i;
	if (upper) {
	    int *nnz = (int *) R_alloc(n, sizeof(int)), nz = Ap[n];
				/* count number of nonzeros in each row */
	    for (i = 0; i < n; i++) nnz[i] = 0;
	    for (i = 0; i < nz; i++) nnz[Ai[i]]++;
	    for (i = 0; i < nf; i++) {
		int j, p2 = Gp[i+1], target = nf - i;
		for (j = Gp[i]; j < p2; j++) {
		    if (nnz[j] != target) return 0;
		}
	    }
	} else {		/* lower triangle - the easy case */
	    for (i = 0; i < nf; i++) {
		int j, p2 = Gp[i+1], target = nf - i;
		for (j = Gp[i]; j < p2; j++) {
		    if ((Ap[j+1] - Ap[j]) != target)
			return 0;
		}
	    }
	}
    }
    return 1;
}

/** 
 * Determine if a fill-reducing permutation is needed for the pairwise
 * cross-tabulation matrix.  If so, determine such a permutation
 * (using Metis) then separate the groups.
 * 
 * @param ctab pointer to a pairwise cross-tabulation object
 * 
 * @return pointer to an integer R vector.
 */

SEXP ctab_permute(SEXP ctab)
{
    SEXP val, GpSl = GET_SLOT(ctab, Matrix_GpSym);
    int *Ai = INTEGER(GET_SLOT(ctab, Matrix_iSym)),
	*Ap = INTEGER(GET_SLOT(ctab, Matrix_pSym)),
	*Gp = INTEGER(GpSl),
	*perm,
	*work,
	i,
	j,
	n = INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],
	nf = length(GpSl) - 1,
	pos;

    if (ctab_isNested(n, nf, 1, Ap, Ai, Gp))
	return allocVector(INTSXP, 0);
    val =  allocVector(INTSXP, n);
    perm = INTEGER(val);
    work = (int *) R_alloc(n, sizeof(int));
    ssc_metis_order(n, Ap, Ai, work, perm);	/* perm gets inverse perm */
    /* work now contains desired permutation but with groups scrambled */

    /* copy work into perm preserving the order of the groups */
    pos = 0;		/* position in new permutation */
    for (i = 0; i < nf; i++) {
	for (j = 0; j < n; j++) {
	    int jj = work[j];
	    if (Gp[i] <= jj && jj < Gp[i+1]) {
		perm[pos] = jj;
		pos++;
	    }
	}
    }
    return val;
}

static
void ssclme_copy_ctab(int nf, const int nc[], SEXP ctab, SEXP ssc)
{
    int *snc, i, copyonly = 1;

    for (i = 0; i < nf; i++) {
	if (nc[i] > 1) copyonly = 0;
    }
    if (copyonly) {
	SET_SLOT(ssc, Matrix_pSym, duplicate(GET_SLOT(ctab, Matrix_pSym)));
	SET_SLOT(ssc, Matrix_iSym, duplicate(GET_SLOT(ctab, Matrix_iSym)));
	SET_SLOT(ssc, Matrix_xSym, duplicate(GET_SLOT(ctab, Matrix_xSym)));
	SET_SLOT(ssc, Matrix_DimSym,
		 duplicate(GET_SLOT(ctab, Matrix_DimSym)));
	SET_SLOT(ssc, Matrix_GpSym, duplicate(GET_SLOT(ctab, Matrix_GpSym)));
    } else {
	int
	    *GpIn = INTEGER(GET_SLOT(ctab, Matrix_GpSym)),
	    *GpOut,
	    *AiIn = INTEGER(GET_SLOT(ctab, Matrix_iSym)),
	    *AiOut,
	    *ApIn = INTEGER(GET_SLOT(ctab, Matrix_pSym)),
	    *ApOut,
	    nIn = GpIn[nf], nOut, nzOut,
	    *dims,
	    *map = Calloc(nIn + 1, int), /* column number map */
	    *ncc = Calloc(nIn, int); /* number of columns out for this
				      * col in */

	SET_SLOT(ssc, Matrix_GpSym, allocVector(INTSXP, nf + 1));
	GpOut = INTEGER(GET_SLOT(ssc, Matrix_GpSym));
	map[0] = GpOut[0] = 0;
	for (i = 0; i < nf; i++) {
	    int j, GpIni = GpIn[i], GpInip1 = GpIn[i+1], nci = nc[i];
	    GpOut[i+1] = GpOut[i] + (GpInip1 - GpIni) * nci;
	    for (j = GpIni; j < GpInip1; j++) {
		ncc[j] = nci;
		map[j+1] = map[j] + nci;
	    }
	}
	nOut = GpOut[nf];	/* size of output matrix */
	SET_SLOT(ssc, Matrix_DimSym,
		 duplicate(GET_SLOT(ctab, Matrix_DimSym)));
	dims = INTEGER(GET_SLOT(ssc, Matrix_DimSym));
	dims[0] = dims[1] = nOut;
	SET_SLOT(ssc, Matrix_pSym, allocVector(INTSXP, nOut + 1));
	ApOut = INTEGER(GET_SLOT(ssc, Matrix_pSym));
	ApOut[0] = 0;
	for (i = 0; i < nf; i++) { /* determine the column pointers */
	    int j, jout = GpOut[i], nci = nc[i], p2 = GpIn[i+1];
	    for (j = GpIn[i]; j < p2; j++) {
		int k, nj = 0, p3 = ApIn[j+1];
		for (k = ApIn[j]; k < p3; k++) {
		    nj += ncc[AiIn[k]];
		}
		nj -= nci - 1;
		ApOut[jout+1] = ApOut[jout] + nj;
		jout++;
		for (k = 1; k < nci; k++) {
		    ApOut[jout+1] = ApOut[jout] + nj + k;
		    jout++;
		}
	    }
	}
	nzOut = ApOut[nOut];	/* number of non-zeros in output */
	SET_SLOT(ssc, Matrix_xSym, allocVector(REALSXP, nzOut));
	memset(REAL(GET_SLOT(ssc, Matrix_xSym)), 0,
	       sizeof(double) * nzOut);
	SET_SLOT(ssc, Matrix_iSym, allocVector(INTSXP, nzOut));
	AiOut = INTEGER(GET_SLOT(ssc, Matrix_iSym));
	for (i = 0; i < nf; i++) { /* fill in the rows */
	    int j, jj, nci = nc[i], p2 = GpIn[i+1];
	    for (j = GpIn[i]; j < p2; j++) { /* first col for input col */
		int ii = AiIn[j], mj = map[j], ncci = ncc[ii],
		    pos = ApOut[mj];
		AiOut[pos++] = map[ii];
		if (ii < j) {	/* above the diagonal */
		    for (jj = 1; jj < ncci; jj++) {
			AiOut[pos+1] = AiOut[pos] + 1;
			pos++;
		    }
		}
		for (jj = 1; jj < nci; jj++) { /* repeat the column adding 
						* another diagonal element */
		    int mjj = mj + jj, pj = ApOut[mjj], pjm1 = ApOut[mjj-1];
		    Memcpy(AiOut + pj, AiOut + pjm1, pj - pjm1);
		    AiOut[ApOut[mjj + 1] - 1] = mjj; /* maybe mjj-1? */
		}
	    }
	}
	Free(map); Free(ncc);
    }
    SET_SLOT(ssc, Matrix_ncSym, allocVector(INTSXP, nf + 2));
    snc = INTEGER(GET_SLOT(ssc, Matrix_ncSym));
    for (i = 0; i <= nf; i++) {
	snc[i] = nc[i];
    }
}

void ssclme_fill_LIp(int n, const int Parent[], int LIp[])
{
    int *sz = Calloc(n, int), i;
    for (i = n - 1; i >= 0; i--) {
	sz[i] = (Parent[i] < 0) ? 0 : 1 + sz[Parent[i]];
    }
    LIp[0] = 0;
    for (i = 0; i < n; i++) LIp[i+1] = LIp[i] + sz[i];
    Free(sz);
}

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]);
	}
    }
}

SEXP
ssclme_create(SEXP facs, SEXP ncv, SEXP threshold)
{
    SEXP ctab, nms, ssc, tmp,
	val = PROTECT(allocVector(VECSXP, 2)),
	dd = PROTECT(allocVector(INTSXP, 3));	/* dimensions of 3-D arrays */
    int *Ai, *Ap, *Gp, *LIp, *Lp, *Parent,
	*nc, Lnz, i, nf = length(facs), nzcol, pp1,
	*dims = INTEGER(dd);

    if (length(ncv) != (nf + 1))
	error("length of nc (%d) should be length of facs (%d) + 1",
	      length(ncv), nf);
    SET_VECTOR_ELT(val, 0, NEW_OBJECT(MAKE_CLASS("ssclme")));
    ssc = VECTOR_ELT(val, 0);
				/* Pairwise cross-tabulation */
    ctab = PROTECT(sscCrosstab(facs, ScalarLogical(1)));
    SET_VECTOR_ELT(val, 1, ctab_permute(ctab));
    if (length(VECTOR_ELT(val, 1)) > 0) {/* Fill-reducing permutation */
	ssc_symbolic_permute(INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],
			     1, INTEGER(VECTOR_ELT(val, 1)),
			     INTEGER(GET_SLOT(ctab, Matrix_pSym)),
			     INTEGER(GET_SLOT(ctab, Matrix_iSym)));
    }
    ssclme_copy_ctab(nf, INTEGER(ncv), ctab, ssc);
    UNPROTECT(1);		/* ctab */

    nzcol = INTEGER(GET_SLOT(ssc, Matrix_DimSym))[1];
    Gp = INTEGER(GET_SLOT(ssc, Matrix_GpSym));
    Ap = INTEGER(GET_SLOT(ssc, Matrix_pSym));
    Ai = INTEGER(GET_SLOT(ssc, Matrix_iSym));
    nc = INTEGER(GET_SLOT(ssc, Matrix_ncSym));
    nc[nf + 1] = length(VECTOR_ELT(facs, 0)); /* number of observations */
				/* Create slots */
    pp1 = nc[nf];
    SET_SLOT(ssc, Matrix_XtXSym, allocMatrix(REALSXP, pp1, pp1));
    SET_SLOT(ssc, Matrix_RXXSym, allocMatrix(REALSXP, pp1, pp1));
    SET_SLOT(ssc, Matrix_ZtXSym, allocMatrix(REALSXP, nzcol, pp1));
    SET_SLOT(ssc, Matrix_RZXSym, allocMatrix(REALSXP, nzcol, pp1));
				/* Zero symmetric matrices (cosmetic) */
    memset(REAL(GET_SLOT(ssc, Matrix_XtXSym)), 0,
	   sizeof(double) * pp1 * pp1); 
    memset(REAL(GET_SLOT(ssc, Matrix_RXXSym)), 0,
	   sizeof(double) * pp1 * pp1);
    SET_SLOT(ssc, Matrix_LpSym, allocVector(INTSXP, nzcol + 1));
    Lp = INTEGER(GET_SLOT(ssc, Matrix_LpSym));
    SET_SLOT(ssc, Matrix_ParentSym, allocVector(INTSXP, nzcol));
    Parent = INTEGER(GET_SLOT(ssc, Matrix_ParentSym));
    SET_SLOT(ssc, Matrix_DSym, allocVector(REALSXP, nzcol));
    SET_SLOT(ssc, Matrix_DIsqrtSym, allocVector(REALSXP, nzcol));
    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 */
    Lnz = Lp[nzcol];
    SET_SLOT(ssc, Matrix_LiSym, allocVector(INTSXP, Lnz));
    SET_SLOT(ssc, Matrix_LxSym, allocVector(REALSXP, Lnz));
    SET_SLOT(ssc, Matrix_OmegaSym, allocVector(VECSXP, nf));
    tmp = GET_SLOT(ssc, Matrix_OmegaSym);
    setAttrib(tmp, R_NamesSymbol, getAttrib(facs, R_NamesSymbol));
    for (i = 0; i < nf; i++) {
	SET_VECTOR_ELT(tmp, i, allocMatrix(REALSXP, nc[i], nc[i]));
	memset(REAL(VECTOR_ELT(tmp, i)), 0,
	       sizeof(double) * nc[i] * nc[i]);
    }
    SET_SLOT(ssc, Matrix_devianceSym, allocVector(REALSXP, 2));
    tmp = GET_SLOT(ssc, Matrix_devianceSym);
    setAttrib(tmp, R_NamesSymbol, allocVector(STRSXP, 2));
    nms = getAttrib(tmp, R_NamesSymbol);
    SET_STRING_ELT(nms, 0, mkChar("ML"));
    SET_STRING_ELT(nms, 1, mkChar("REML"));
    SET_SLOT(ssc, Matrix_devCompSym, allocVector(REALSXP, 4));
    SET_SLOT(ssc, Matrix_statusSym, allocVector(LGLSXP, 2));
    tmp = GET_SLOT(ssc, Matrix_statusSym);
    LOGICAL(tmp)[0] = LOGICAL(tmp)[1] = 0;
    setAttrib(tmp, R_NamesSymbol, allocVector(STRSXP, 2));
    nms = getAttrib(tmp, R_NamesSymbol);
    SET_STRING_ELT(nms, 0, mkChar("factored"));
    SET_STRING_ELT(nms, 1, mkChar("inverted"));
    SET_SLOT(ssc, Matrix_bVarSym, allocVector(VECSXP, nf));
    tmp = GET_SLOT(ssc, Matrix_bVarSym);
    setAttrib(tmp, R_NamesSymbol, getAttrib(facs, R_NamesSymbol));
    for (i = 0; i < nf; i++) {
	int nci = nc[i], mi = (Gp[i+1] - Gp[i])/nc[i];

	dims[0] = dims[1] = nci;
	dims[2] = mi;
	SET_VECTOR_ELT(tmp, i, allocArray(REALSXP, dd));
	memset(REAL(VECTOR_ELT(tmp, i)), 0,
	       sizeof(double) * nci * nci * mi);
    }
    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);
    }
    UNPROTECT(2);
    return val;
}

static
void bVj_to_A(int ncj, int Gpj, int Gpjp, const double bVj[],
	      const int Ap[], const int Ai[], double Ax[])
{
    int i, diag, k;
    for (i = Gpj; i < Gpjp; i += ncj) {
	for (k = 0; k < ncj; k++) {
	    diag = Ap[i + k + 1] - 1;
	    if (Ai[diag] != i+k)
		error("Expected Ai[%d] to be %d (i.e on diagonal) not %d",
		      diag, i+k, Ai[diag]);
	    Memcpy(Ax + diag - k, bVj + (i+k-Gpj)*ncj, k + 1);
	}
    }
}

SEXP
ssclme_transfer_dimnames(SEXP x, SEXP facs, SEXP mmats)
{
    SEXP bVar = GET_SLOT(x, Matrix_bVarSym),
	nms2 = PROTECT(allocVector(VECSXP, 2)),
	nms3 = PROTECT(allocVector(VECSXP, 3));
    int i, nf = length(mmats) - 1;
    SEXP xcols = VECTOR_ELT(GetArrayDimnames(VECTOR_ELT(mmats, nf)), 1);

    for (i = 0; i < nf; i++) {
	SEXP cnms = VECTOR_ELT(GetArrayDimnames(VECTOR_ELT(mmats, i)), 1);
	SET_VECTOR_ELT(nms3, 0, cnms);
	SET_VECTOR_ELT(nms3, 1, cnms);
	SET_VECTOR_ELT(nms3, 2,
		       getAttrib(VECTOR_ELT(facs, i), R_LevelsSymbol));
	dimnamesgets(VECTOR_ELT(bVar, i), duplicate(nms3));
    }
    SET_VECTOR_ELT(nms2, 0, xcols);
    SET_VECTOR_ELT(nms2, 1, xcols);
    dimnamesgets(GET_SLOT(x, Matrix_XtXSym), nms2);
    dimnamesgets(GET_SLOT(x, Matrix_RXXSym), nms2);
    UNPROTECT(2);
    return R_NilValue;
}

SEXP
ssclme_update_mm(SEXP x, SEXP facs, SEXP mmats)
{
    SEXP bVar = GET_SLOT(x, Matrix_bVarSym);
    int
	*Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
	*Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
	*Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
	*nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
	i, j, k,
	ione = 1,
	nf = length(mmats) - 1,
	nobs = nc[nf + 1],
	nzcol = Gp[nf],
	nz = Ap[nzcol],
	pp1 = nc[nf];
    double
	**Z = Calloc(nf + 1, double *),
	*Ax = REAL(GET_SLOT(x, Matrix_xSym)),
	*XtX = REAL(GET_SLOT(x, Matrix_XtXSym)),
	*ZtX = REAL(GET_SLOT(x, Matrix_ZtXSym)),
	one = 1.0,
	zero = 0.0;

    for (i = 0; i <= nf; i++) {
	int *dims = INTEGER(getAttrib(VECTOR_ELT(mmats, i), R_DimSymbol)),
	    nci = nc[i];
	if (nobs != dims[0])
	    error("Expected %d rows in the %d'th model matrix. Got %d",
		  nobs, i+1, dims[0]);
	if (nci != dims[1])
	    error("Expected %d columns in the %d'th model matrix. Got %d",
		  nci, i+1, dims[1]);
	Z[i] = REAL(VECTOR_ELT(mmats, i));
    }
				/* Create XtX - X is Z[nf] */
    F77_CALL(dsyrk)("U", "T", nc+nf, &nobs, &one,
		    Z[nf], &nobs, &zero, XtX, nc + nf);
				/* Zero the accumulators */
    memset((void *) ZtX, 0, sizeof(double) * nzcol * pp1);
    memset((void *) Ax, 0, sizeof(double) * nz);
    for (j = 0; j < nf; j++) { /* Create ZtX */
	int *fpj = INTEGER(VECTOR_ELT(facs, j)), ncj = nc[j],
	    Ncj = ncj > 1;
	double
	    *bVj = REAL(VECTOR_ELT(bVar, j)),
	    *Zj = Z[j],
	    *zxj = ZtX + Gp[j];

	if (Ncj) {		/* bVj will accumulate Z'Z blocks */
	    memset(bVj, 0, sizeof(double) * ncj * (Gp[j+1]-Gp[j]));
	}
	for (i = 0; i < nobs; i++) { /* accumulate diagonal of ZtZ */
	    int fpji = fpj[i] - 1, /* factor indices are 1-based */
		dind = Ap[Gp[j] + fpji * ncj + 1] - 1;
	    if (Ai[dind] != (Gp[j] + fpji * ncj))
		error("logic error in ssclme_update_mm");
	    if (Ncj) {		/* use bVar to accumulate */
		F77_CALL(dsyrk)("U", "T", &ncj, &ione, &one, Zj+i,
				&nobs, &one, bVj + fpji*ncj*ncj, &ncj);
	    } else {		/* update scalars directly */
		Ax[dind] += Zj[i] * Zj[i];
	    }
				/* update rows of Z'X */
	    F77_CALL(dgemm)("T", "N", &ncj, &pp1, &ione, &one,
			    Zj + i, &nobs, Z[nf] + i, &nobs,
			    &one, zxj + fpji * ncj, &nzcol);
	}
	if (Ncj) bVj_to_A(ncj, Gp[j], Gp[j+1], bVj, Ap, Ai, Ax);
	for (k = j+1; k < nf; k++) { /* off-diagonals */
	    int *fpk = INTEGER(VECTOR_ELT(facs, k)),
		*Apk = Ap + Gp[k],
		nck = nc[k];
	    double
		*Zk = Z[k];

	    for (i = 0; i < nobs; i++) {
		int ii, ind = -1, fpji = fpj[i] - 1,
		    row = Gp[j] + fpji * ncj,
		    fpki = fpk[i] - 1,
		    lastind = Apk[fpki + 1];
		for (ii = Apk[fpki]; ii < lastind; ii++) {
		    if (Ai[ii] == row) {
			ind = ii;
			break;
		    }
		}
		if (ind < 0) error("logic error in ssclme_update_mm");
		if (Ncj || nck > 1) {
				/* FIXME: run a loop to update */
		    error("code not yet written");
		} else {	/* update scalars directly */
		    Ax[ind] += Zj[fpji] * Zk[fpki];
		}
	    }
	}
    }
    Free(Z);
    ssclme_transfer_dimnames(x, facs, mmats);
    return R_NilValue;
}

SEXP ssclme_inflate_and_factor(SEXP x)
{
    SEXP
	GpSlot = GET_SLOT(x, Matrix_GpSym),
	Omega = GET_SLOT(x, Matrix_OmegaSym);
    int n = INTEGER(GET_SLOT(x, Matrix_DimSym))[1];
    int
	*Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
	*Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
	*Flag = Calloc(n, int),
	*Gp = INTEGER(GpSlot),
	*Lnz = Calloc(n, int),
	*Pattern = Calloc(n, int),
	*nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
	j,
	nf = length(GpSlot) - 1;
    double
	*D = REAL(GET_SLOT(x, Matrix_DSym)),
	*DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
	*Y = Calloc(n, double),
	*xcp = Calloc(Ap[n], double);

    Memcpy(xcp, REAL(GET_SLOT(x, Matrix_xSym)), Ap[n]);
    for (j = 0; j < nf; j++) {
	int  diag, i, ii, k, G2 = Gp[j + 1], ncj = nc[j];
	double *omgj = REAL(VECTOR_ELT(Omega, j));
	
	for (i = Gp[j]; i < G2; i += ncj) {
	    for (k = 0; k < ncj; k++) {
		diag = Ap[i + k + 1] - 1;
		if (Ai[diag] != i+k)
		    error("Expected Ai[%d] to be %d (i.e on diagonal) not %d",
			  diag, i+k, Ai[diag]);
		for (ii = 0; ii <= k; ii++) {
		    xcp[diag + ii - k] += omgj[k*ncj + ii];
		}
	    }
	}
    }
    j = ldl_numeric(n, Ap, Ai, xcp,
		    INTEGER(GET_SLOT(x, Matrix_LpSym)),
		    INTEGER(GET_SLOT(x, Matrix_ParentSym)),
		    Lnz, INTEGER(GET_SLOT(x, Matrix_LiSym)),
		    REAL(GET_SLOT(x, Matrix_LxSym)),
		    D, Y, Pattern, Flag,
		    (int *) NULL, (int *) NULL); /* P & Pinv */
    if (j != n)
	error("rank deficiency of ZtZ+W detected at column %d",
	      j + 1);
    for (j = 0; j < n; j++) DIsqrt[j] = 1./sqrt(D[j]);
    Free(Lnz); Free(Flag); Free(Pattern); Free(Y); Free(xcp);
    return R_NilValue;
}

SEXP ssclme_factor(SEXP x)
{
    int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
    
    if (!status[0]) {
	SEXP
	    GpSlot = GET_SLOT(x, Matrix_GpSym),
	    Omega = GET_SLOT(x, Matrix_OmegaSym);
	int
	    *Gp = INTEGER(GpSlot),
	    *Li = INTEGER(GET_SLOT(x, Matrix_LiSym)),
	    *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
	    *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
	    i,
	    n = INTEGER(GET_SLOT(x, Matrix_DimSym))[1],
	    nf = length(GpSlot) - 1,
	    nobs = nc[nf + 1],
	    nreml = nobs + 1 - nc[nf],
	    pp1 = nc[nf],
	    pp2 = pp1 + 1;
	double
	    *D = REAL(GET_SLOT(x, Matrix_DSym)),
	    *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
	    *Lx = REAL(GET_SLOT(x, Matrix_LxSym)),
	    *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
	    *RZX = REAL(GET_SLOT(x, Matrix_RZXSym)),
	    *dcmp = REAL(getAttrib(x, Matrix_devCompSym)),
	    *deviance = REAL(getAttrib(x, Matrix_devianceSym)),
	    minus1 = -1.,
	    one = 1.;
	
	ssclme_inflate_and_factor(x);
				/* Accumulate logdet of ZtZ+W */
	dcmp[0] = dcmp[1] = dcmp[2] = dcmp[3] = 0.;
	for (i = 0; i < n; i++) dcmp[0] += log(D[i]);
				/* Accumulate logdet of W */
	for (i = 0; i < nf; i++) {
	    int nci = nc[i],
		mi = (Gp[i+1] - Gp[i])/nci;
	
	    if (nci < 2) {
		dcmp[1] += mi * log(REAL(VECTOR_ELT(Omega, i))[0]);
	    } else {
		int j;
		double
		    *tmp = Calloc(nci * nci, double),
		    accum = 0.;
		F77_CALL(dpotrf)("U", &nci,
				 Memcpy(tmp, REAL(VECTOR_ELT(Omega, i)),
					nci * nci),
				 &nci, &j);
		if (j) 
		    error("Omega[%d] is not positive definite", i + 1);
		for (j = 0; j < nci; j++) {
		    accum += 2 * log(tmp[j * (nci + 1)]);
		}
		dcmp[1] += mi * accum;
		Free(tmp);
	    }
	}
				/* ldl_lsolve on Z'X */
	Memcpy(RZX, REAL(GET_SLOT(x, Matrix_ZtXSym)), n * pp1);
	for (i = 0; i < pp1; i++) {
	    int j;
	    double *RZXi = RZX + i * n;
	    ldl_lsolve(n, RZXi, Lp, Li, Lx);
	    for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
	}
				/* downdate and factor X'X */
	Memcpy(RXX, REAL(GET_SLOT(x, Matrix_XtXSym)), pp1 * pp1);
	F77_CALL(dsyrk)("U", "T", &pp1, &n, &minus1,
			RZX, &n, &one, RXX, &pp1);
	F77_CALL(dpotrf)("U", &pp1, RXX, &pp1, &i);
	if (i)
	    error("DPOTRF returned error code %d", i);
				/* logdet of RXX */
	for (i = 0; i < (pp1 - 1); i++)
	    dcmp[2] += 2 * log(RXX[i*pp2]);
				/* logdet of Ryy */
	dcmp[3] = 2. * log(RXX[pp1 * pp1 - 1]);
	deviance[0] =		/* ML criterion */
	    dcmp[0] - dcmp[1] + nobs*(1+dcmp[3]+log(2*PI/nobs));
	deviance[1] = dcmp[0] - dcmp[1] + /* REML */
	    dcmp[2] + nreml * (1. + dcmp[3] + log(2. * PI/nreml));
	status[0] = 1;
	status[1] = 0;
    }
    return R_NilValue;
}

static
int ldl_update_ind(int probe, int start, const int ind[])
{
    while (ind[start] < probe) start++;
    if (ind[start] > probe) error("logic error in ldl_inverse");
    return start;
}

/** 
 * Create the inverse of L and update the diagonal blocks of the inverse
 * of LDL' (=Z'Z+W)
 * 
 * @param x pointer to an ssclme object
 *
 * @return R_NilValue (x is updated in place)

 */
static
SEXP ldl_inverse(SEXP x)
{
    SEXP
	Gpsl = GET_SLOT(x, Matrix_GpSym),
	LIisl = GET_SLOT(x, Matrix_LIiSym),
	LIpsl = GET_SLOT(x, Matrix_LIpSym),
	bVar = GET_SLOT(x, Matrix_bVarSym);
    int *Gp = INTEGER(Gpsl),
	*Li,
	*LIp = INTEGER(LIpsl), *Lp,
	i,
	nf = length(Gpsl) - 1,
	nzc = length(LIpsl) - 1;
    double
	*DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
	*Lx;
	
    ssclme_factor(x);
    if (LIp[nzc] == 0) {	/* L and LI are the identity */
	for (i = 0; i < nf; i++) {
	    Memcpy(REAL(VECTOR_ELT(bVar, i)), DIsqrt + Gp[i],
		   Gp[i+1] - Gp[i]);
	}
	return R_NilValue;
    }
    Lp = INTEGER(GET_SLOT(x, Matrix_LpSym));
    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;

	double *LIx = REAL(GET_SLOT(x, Matrix_LIxSym)),
	    one = 1., zero = 0.;

	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 (cosmetic) */
	    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)("U", "T", &nci, &rr, &one, tmp, &nr,
				&zero, bVi + (j - G1) * nci, &nci);
		F77_CALL(dpotrf)("U", &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);
	    }
	}
	return R_NilValue;
    }
    if (length(LIisl)) error("logic error in ssclme_ldl_inverse");
    else {		    /* LIi and LIx are too big and not used */
	int *counts = Calloc(nzc, int), info, maxod = -1;
	int *Parent = INTEGER(GET_SLOT(x, Matrix_ParentSym));
	int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym));
	double one = 1.0, zero = 0.;
				/* determine maximum # of off-diagonals */
	for (i = nzc - 1; i >= 0; i--) { /* in a column of L^{-1} */
	    counts[i] = (Parent[i] < 0) ? 0 : 1 + counts[Parent[i]];
	    if (counts[i] > maxod) maxod = counts[i];
	}
	Free(counts);

	for (i = 0; i < nf; i++) {
	    int j, jj, k, kk, nci = nc[i], nr, p, p2, pp,
		m = maxod + nci,
		*ind = Calloc(m, int);
	    double
		*tmp = Calloc(m * nci, double),
		*mpt = REAL(VECTOR_ELT(bVar, i));
	    
	    for (j = Gp[i]; j < Gp[i+1]; j += nci) {
		memset((void *) tmp, 0, sizeof(double) * m * nci);
		
		kk = 0;		/* ind holds indices of non-zeros */
		jj = j;		/* in this block of columns */
		while (jj >= 0) {
		    ind[kk++] = jj;
		    jj = Parent[jj];
		}
		nr = kk;       /* number of non-zeros in this block */
		while (kk < m) ind[kk++] = nzc; /* placeholders */

		for (k = 0; k < nci; k++) {
		    double *ccol = tmp + k * nr;
		    
		    ccol[k] = 1.;
		    kk = k;
		    for (jj = j + k; jj >= 0; jj = Parent[jj]) {
			p2 = Lp[jj+1];
			pp = kk;
			for (p = Lp[jj]; p < p2; p++) {
			    pp = ldl_update_ind(Li[p], pp, ind);
			    ccol[pp] -= Lx[p] * ccol[kk];
			}
		    }
		}
				/* scale rows */
		for (kk = 0; kk < nr; kk++) {
		    for (k = 0; k < nci; k++) {
			tmp[k * nr + kk] *= DIsqrt[ind[kk]];
		    }
		}
		F77_CALL(dsyrk)("U", "T", &nci, &nr, &one, tmp, &nr,
				&zero, mpt + (j - Gp[i])*nci, &nci);
		F77_CALL(dpotrf)("U", &nci, mpt + (j - Gp[i])*nci,
				 &nci, &info);
		if (info)	/* should never happen */
		    error(
			"Rank deficient variance matrix at group %d, level %d",
			i + 1, j + 1);
	    }
	    Free(tmp); Free(ind);
	}
    }
    return R_NilValue;
}

SEXP ssclme_invert(SEXP x)
{
    int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
    if (!status[0]) ssclme_factor(x);
    if (!status[1]) {
	SEXP
	    RZXsl = GET_SLOT(x, Matrix_RZXSym);
	int
	    *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
	    *Li = INTEGER(GET_SLOT(x, Matrix_LiSym)),
	    *Lp = INTEGER(GET_SLOT(x, Matrix_LpSym)),
	    i,
	    n = dims[0],
	    pp1 = dims[1];
	double
	    *DIsqrt = REAL(GET_SLOT(x, Matrix_DIsqrtSym)),
	    *Lx = REAL(GET_SLOT(x, Matrix_LxSym)),
	    *RXX = REAL(GET_SLOT(x, Matrix_RXXSym)),
	    *RZX = REAL(RZXsl),
	    one = 1.;

	F77_CALL(dtrtri)("U", "N", &pp1, RXX, &pp1, &i);
	if (i)
	    error("DTRTRI returned error code %d", i);
	F77_CALL(dtrmm)("R", "U", "N", "N", &n, &pp1, &one,
			RXX, &pp1, RZX, &n);
	for (i = 0; i < pp1; i++) {
	    int j; double *RZXi = RZX + i * n;
	    for (j = 0; j < n; j++) RZXi[j] *= DIsqrt[j];
	    ldl_ltsolve(n, RZXi, Lp, Li, Lx);
	}
	ldl_inverse(x);
	status[1] = 1;
    }
    return R_NilValue;
}

SEXP ssclme_initial(SEXP x)
{
    SEXP Gpsl = GET_SLOT(x, Matrix_GpSym),
	Omg = GET_SLOT(x, Matrix_OmegaSym);
    int *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
	*Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
	*Gp = INTEGER(Gpsl),
	*nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
	*status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
	i, nf = length(Gpsl) - 1;
	
    double *Ax = REAL(GET_SLOT(x, Matrix_xSym));

    for (i = 0; i < nf; i++) {
	int
	    Gpi = Gp[i],
	    j, k,
	    nci = nc[i],
	    ncip1 = nci + 1,
	    p2 = Gp[i+1];
	double
	    mi = 0.375 * ((double) nci)/((double) (p2 - Gpi)),
	    *mm = REAL(VECTOR_ELT(Omg, i));

	memset((void *) mm, 0, sizeof(double) * nci * nci);
	for (j = Gpi; j < p2; j += nci) {
	    for (k = 0; k < nci; k++) {
		int jk = j+k, jj = Ap[jk+1] - 1;
		if (Ai[jj] != jk) error("malformed ZtZ structure");
		mm[k * ncip1] += Ax[jj] * mi;
	    }
	}
    }
    status[0] = status[1] = 0;
    return R_NilValue;
}

/** 
 * Extract the conditional estimates of the fixed effects
 * 
 * @param x Pointer to an ssclme object
 * 
 * @return a numeric vector containing the conditional estimates of
 * the fixed effects
 */
SEXP ssclme_fixef(SEXP x)
{
    SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);
    int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1];
    int j, p = pp1 - 1;
    SEXP val = PROTECT(allocVector(REALSXP, p));
    double
	*RXX = REAL(RXXsl),
	*beta = REAL(val),
	nryyinv;		/* negative ryy-inverse */

    ssclme_invert(x);
    Memcpy(beta, RXX + p * pp1, p);
    nryyinv = -RXX[pp1*pp1 - 1];
    for (j = 0; j < p; j++) beta[j] /= nryyinv;
    UNPROTECT(1);
    return val;
}

/** 
 * Extract the conditional modes of the random effects.
 * 
 * @param x Pointer to an ssclme object
 * 
 * @return a vector containing the conditional modes of the random effects
 */
SEXP ssclme_ranef(SEXP x)
{
    SEXP RZXsl = GET_SLOT(x, Matrix_RZXSym),
	GpSl = GET_SLOT(x, Matrix_GpSym);
    int *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
	*Gp = INTEGER(GpSl),
	*nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
	i, j,
	n = dims[0],
	nf = length(GpSl) - 1,
	pp1 = dims[1],
	p = pp1 - 1;
    SEXP val = PROTECT(allocVector(VECSXP, nf));
    double
	*b = REAL(RZXsl) + n * p,
	ryyinv;		/* ryy-inverse */

    ssclme_invert(x);
    ryyinv = REAL(GET_SLOT(x, Matrix_RXXSym))[pp1*pp1 - 1];
    for (i = 0; i < nf; i++) {
	int nci = nc[i], Mi = Gp[i+1] - Gp[i];
	double *mm;
	
	SET_VECTOR_ELT(val, i, allocMatrix(REALSXP, nci, Mi/nci));
	mm = Memcpy(REAL(VECTOR_ELT(val, i)), b, Mi);
	b += Mi;
	for (j = 0; j < Mi; j++) mm[j] /= ryyinv;
    }
    UNPROTECT(1);
    return val;
}

/** 
 * Extract the ML or REML conditional estimate of sigma
 * 
 * @param x pointer to an ssclme object
 * @param REML logical scalar - TRUE if REML estimates are requested
 * 
 * @return numeric scalar 
 */
SEXP ssclme_sigma(SEXP x, SEXP REML)
{
    SEXP RXXsl = GET_SLOT(x, Matrix_RXXSym);
    int pp1 = INTEGER(getAttrib(RXXsl, R_DimSymbol))[1],
	nobs = INTEGER(GET_SLOT(x, Matrix_ncSym))
	[length(GET_SLOT(x, Matrix_GpSym))];
   
    ssclme_invert(x);
    return ScalarReal(1./(REAL(RXXsl)[pp1*pp1 - 1] *
			  sqrt((double)(asLogical(REML) ?
					nobs + 1 - pp1 : nobs))));
}

static
int coef_length(int nf, const int nc[])
{
    int i, ans = 0;
    for (i = 0; i < nf; i++) ans += (nc[i] * (nc[i] + 1))/2;
    return ans;
}

/** 
 * Extract the upper triangles of the Omega matrices.
 * (These aren't "coefficients" but the extractor is
 * called coef for historical reasons.)
 * 
 * @param x pointer to an ssclme object
 * 
 * @return numeric vector of the values in the upper triangles of the
 * Omega matrices
 */
SEXP ssclme_coef(SEXP x)
{
    SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
    int	*nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
	i, nf = length(Omega), vind;
    SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
    double *vv = REAL(val);

    vind = 0;
    for (i = 0; i < nf; i++) {
	int nci = nc[i];
	if (nci == 1) {
	    vv[vind++] = REAL(VECTOR_ELT(Omega, i))[0];
	} else {
	    int j, k, odind = vind + nci, ncip1 = nci + 1;
	    double *omgi = REAL(VECTOR_ELT(Omega, i));
	    
	    for (j = 0; j < nci; j++) {
		vv[vind++] = omgi[j * ncip1];
		for (k = j + 1; k < nci; k++) {
		    vv[odind++] = omgi[k*nci + j];
		}
	    }
	    vind = odind;
	}
    }
    UNPROTECT(1);
    return val;
}

/** 
 * Extract the unconstrained parameters that determine the
 * Omega matrices. (Called coef for historical reasons.)
 * 
 * @param x pointer to an ssclme object
 * 
 * @return numeric vector of unconstrained parameters that determine the
 * Omega matrices
 */
SEXP ssclme_coefUnc(SEXP x)
{
    SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
    int	*nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
	i, nf = length(Omega), vind;
    SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
    double *vv = REAL(val);

    vind = 0;
    for (i = 0; i < nf; i++) {
	int  nci = nc[i];
	if (nci == 1) {
	    vv[vind++] = log(REAL(VECTOR_ELT(Omega, i))[0]);
	} else {
	    int j, k, ncip1 = nci + 1, ncisq = nci * nci;
	    double *tmp = Memcpy(Calloc(ncisq, double), 
				 REAL(VECTOR_ELT(Omega, i)), ncisq);
	    F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
	    if (j)		/* should never happen */
		error("DPOTRF returned error code %d on Omega[[%d]]",
		      j, i+1);
	    for (j = 0; j < nci; j++) {
		double diagj = tmp[j * ncip1];
		vv[vind++] = 2. * log(diagj);
		for (k = j + 1; k < nci; k++) {
		    tmp[j + k * nci] /= diagj;
		}
	    }
	    for (j = 0; j < nci; j++) {
		for (k = j + 1; k < nci; k++) {
		    vv[vind++] = tmp[j + k * nci];
		}
	    }
	    Free(tmp);
	}
    }
    UNPROTECT(1);
    return val;
}

/** 
 * Assign the Omega matrices from the unconstrained parameterization.
 * 
 * @param x pointer to an ssclme object
 * @param coef pointer to an numeric vector of appropriate length
 * 
 * @return R_NilValue
 */
SEXP ssclme_coefGetsUnc(SEXP x, SEXP coef)
{
    SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
    int	*nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
	cind, i, nf = length(Omega),
	*status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
    double *cc = REAL(coef);

    if (length(coef) != coef_length(nf, nc) || !isReal(coef))
	error("coef must be a numeric vector of length %d",
	      coef_length(nf, nc));
    cind = 0;
    for (i = 0; i < nf; i++) {
	int nci = nc[i];
	if (nci == 1) {
	    REAL(VECTOR_ELT(Omega, i))[0] = exp(cc[cind++]);
	} else {
	    int odind = cind + nci, /* off-diagonal index */
		j, k,
		ncip1 = nci + 1,
		ncisq = nci * nci;
	    double
		*omgi = REAL(VECTOR_ELT(Omega, i)),
		*tmp = Calloc(ncisq, double),
		diagj, one = 1., zero = 0.;

	    memset(omgi, 0, sizeof(double) * ncisq);
	    for (j = 0; j < nci; j++) {
		tmp[j * ncip1] = diagj = exp(cc[cind++]/2.);
		for (k = j + 1; k < nci; k++) {
		    tmp[k*nci + j] = cc[odind++] * diagj;
		}
	    }
	    F77_CALL(dsyrk)("U", "T", &nci, &nci, &one,
			    tmp, &nci, &zero, omgi, &nci);
	    Free(tmp);
	    cind = odind;
	}
    }
    status[0] = status[1] = 0;
    return x;
}

/** 
 * Assign the upper triangles of the Omega matrices.
 * (Called coef for historical reasons.)
 * 
 * @param x pointer to an ssclme object
 * @param coef pointer to an numeric vector of appropriate length
 * 
 * @return R_NilValue
 */
SEXP ssclme_coefGets(SEXP x, SEXP coef)
{
    SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
    int	*nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
	cind, i, nf = length(Omega),
	*status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
    double *cc = REAL(coef);

    if (length(coef) != coef_length(nf, nc) || !isReal(coef))
	error("coef must be a numeric vector of length %d",
	      coef_length(nf, nc));
    cind = 0;
    for (i = 0; i < nf; i++) {
	int nci = nc[i];
	if (nci == 1) {
	    REAL(VECTOR_ELT(Omega, i))[0] = cc[cind++];
	} else {
	    int j, k, odind = cind + nci, ncip1 = nci + 1;
	    double *omgi = REAL(VECTOR_ELT(Omega, i));
	
	    for (j = 0; j < nci; j++) {
		omgi[j * ncip1] = cc[cind++];
		for (k = j + 1; k < nci; k++) {
		    omgi[k*nci + j] = cc[odind++];
		}
	    }
	    cind = odind;
	}
    }
    status[0] = status[1] = 0;
    return x;
}

SEXP ssclme_EMsteps(SEXP x, SEXP nsteps, SEXP REMLp, SEXP verb)
{
    SEXP
	Omega = GET_SLOT(x, Matrix_OmegaSym),
	RZXsl = GET_SLOT(x, Matrix_RZXSym),
	ncsl = GET_SLOT(x, Matrix_ncSym),
	bVar = GET_SLOT(x, Matrix_bVarSym);
    int
	*Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
	*dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
	*nc = INTEGER(ncsl),
	*status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
	REML = asLogical(REMLp),
	i, info, iter,
	n = dims[0],
	nEM = asInteger(nsteps),
	nf = length(ncsl) - 2,
	nobs = nc[nf + 1],
	p,
	pp1 = dims[1],
	verbose = asLogical(verb);
    double
	*RZX = REAL(RZXsl),
	*dev = REAL(GET_SLOT(x, Matrix_devianceSym)),
	*b,
        alpha,
	one = 1.,
	zero = 0.;

    p = pp1 - 1;
    b = RZX + p * n;
    if (verbose) Rprintf("  EM iterations\n");
    for (iter = 0; iter <= nEM; iter++) {
	ssclme_invert(x);
	if (verbose) {
	    SEXP coef = PROTECT(ssclme_coef(x));
	    int lc = length(coef); double *cc = REAL(coef);
	    Rprintf("%3d %.3f", iter, dev[REML ? 1 : 0]);
	    for (i = 0; i < lc; i++) Rprintf(" %#8g", cc[i]);
	    Rprintf("\n");
	    UNPROTECT(1);
	}
	for (i = 0; i < nf; i++) {
	    int ki = Gp[i+1] - Gp[i],
		nci = nc[i],
		mi = ki/nci;
	    double
	        *vali = REAL(VECTOR_ELT(Omega, i));
	    
	    alpha = ((double)(REML?(nobs-p):nobs))/((double)mi);
	    F77_CALL(dsyrk)("U", "N", &nci, &mi,
			    &alpha, b + Gp[i], &nci,
			    &zero, vali, &nci);
	    alpha = 1./((double) mi);
	    F77_CALL(dsyrk)("U", "N", &nci, &ki,
			    &alpha, REAL(VECTOR_ELT(bVar, i)), &nci,
			    &one, vali, &nci);
	    if (REML) {
		int j;
		for (j = 0; j < p; j++) { 
		    F77_CALL(dsyrk)("U", "N", &nci, &mi,
				&alpha, RZX + Gp[i] + j*n, &nci,
				&one, vali, &nci);
		}
	    }
	    F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);
	    if (info) error("DPOTRF returned error code %d", info);
	    F77_CALL(dpotri)("U", &nci, vali, &nci, &info);
	    if (info) error("DPOTRI returned error code %d", info);
	}
	status[0] = status[1] = 0;
    }
    ssclme_factor(x);
    return R_NilValue;
}

SEXP ssclme_gradient(SEXP x, SEXP REMLp, SEXP Uncp)
{
    SEXP
	Omega = GET_SLOT(x, Matrix_OmegaSym),
	RZXsl = GET_SLOT(x, Matrix_RZXSym),
	ncsl = GET_SLOT(x, Matrix_ncSym),
	bVar = GET_SLOT(x, Matrix_bVarSym);
    int
	*Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
	*dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
	*nc = INTEGER(ncsl),
	REML = asLogical(REMLp),
	cind, i, n = dims[0],
	nf = length(Omega),
	nobs, odind, p, pp1 = dims[1],
	uncst = asLogical(Uncp);
    double
	*RZX = REAL(RZXsl),
	*b,
        alpha,
	one = 1.;
    SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));

    nobs = nc[nf + 1];
    p = pp1 - 1;
    b = RZX + p * n;
    ssclme_invert(x);
    cind = 0;
    for (i = 0; i < nf; i++) {
	int j, ki = Gp[i+1] - Gp[i],
	    nci = nc[i], ncip1 = nci + 1, ncisq = nci * nci,
	    mi = ki/nci;
	double
	    *chol = Memcpy(Calloc(ncisq, double),
			   REAL(VECTOR_ELT(Omega, i)), ncisq),
	    *tmp = Calloc(ncisq, double);
	
	    
	F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
	if (j)
	    error("DPOTRF gave error code %d on Omega[[%d]]", j, i + 1);
	Memcpy(tmp, chol, ncisq);
	F77_CALL(dpotri)("U", &nci, tmp, &nci, &j);
	if (j)
	    error("DPOTRI gave error code %d on Omega[[%d]]", j, i + 1);
	alpha = (double) -mi;
	F77_CALL(dsyrk)("U", "N", &nci, &ki,
			&one, REAL(VECTOR_ELT(bVar, i)), &nci,
			&alpha, tmp, &nci);
	alpha = ((double)(REML?(nobs-p):nobs));
	F77_CALL(dsyrk)("U", "N", &nci, &mi,
			&alpha, b + Gp[i], &nci,
			&one, tmp, &nci);
	if (REML) {
	    for (j = 0; j < p; j++) { 
		F77_CALL(dsyrk)("U", "N", &nci, &mi,
				&one, RZX + Gp[i] + j*n, &nci,
				&one, tmp, &nci);
	    }
	}
	if (nci == 1) {
	    REAL(ans)[cind++] = *tmp *
		(uncst ? *REAL(VECTOR_ELT(Omega, i)) : 1.);
	} else {
	    int odind = cind + nci;
	    if (uncst) {
		int kk;
		nlme_symmetrize(tmp, nci);
		for (j = 0; j < nci; j++, cind++) {
		    double *dder = REAL(ans) + cind;
		    *dder = 0.;
		    for (k = j; k < nci; k++) {
			for (kk = j; kk < nci; kk++) {
			    *dder += chol[k*nci + j] * chol[kk*nci + j] *
				tmp[kk * nci + k];
			}
		    }
		    for (k = j + 1; k < nci; k++) {
			REAL(ans)[odind++] = -1.;
		    }			
		}
	    } else {
		for (j = 0; j < nci; j++) {
		    REAL(ans)[cind++] = tmp[j * ncip1];
		    for (k = j + 1; k < nci; k++) {
			REAL(ans)[odind++] = tmp[k*nci + j] * 2.;
		    }
		}
	    }
	    cind = odind;
	}
	Free(tmp); Free(chol);
    }
    UNPROTECT(1);
    return ans;
}

SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats)
{
    SEXP val, b;
    int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
	i, ione = 1, nf = length(facs), nobs, p;
    double *vv, one = 1.0, zero = 0.0;

    if (nf < 1)
	error("null factor list passed to ssclme_fitted");
    nobs = length(VECTOR_ELT(facs, 0));
    val = PROTECT(allocVector(REALSXP, nobs));
    vv = REAL(val);
    p = nc[nf] - 1;
    if (p > 0) {
	F77_CALL(dgemm)("N", "N", &nobs, &ione, &p, &one,
			REAL(VECTOR_ELT(mmats, nf)), &nobs,
			REAL(PROTECT(ssclme_fixef(x))), &p,
			&zero, vv, &nobs);
	UNPROTECT(1);
    } else {
	memset(vv, 0, sizeof(double) * nobs);
    }
    b = PROTECT(ssclme_ranef(x));
    for (i = 0; i < nf; i++) {
	int *ff = INTEGER(VECTOR_ELT(facs, i)), j, nci = nc[i];
	double *bb = REAL(VECTOR_ELT(b, i)),
	    *mm = REAL(VECTOR_ELT(mmats, i));
	for (j = 0; j < nobs; ) {
	    int nn = 1, lev = ff[j];
	    /* check for adjacent rows with same factor level */
	    while (ff[j + nn] == lev) nn++; 
	    F77_CALL(dgemm)("N", "N", &nn, &ione, &nci,
			    &one, mm + j, &nobs,
			    bb + (lev - 1) * nci, &nci,
			    &one, vv + j, &nobs);
	    j += nn;
	}
    }
    UNPROTECT(2);
    return val;
}

SEXP ssclme_variances(SEXP x, SEXP REML)
{
    SEXP Omg = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym))),
	val = PROTECT(allocVector(VECSXP, 2));
    int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
	i, nf = length(Omg);
    double sigmasq;
    

    SET_VECTOR_ELT(val, 0, Omg);
    SET_VECTOR_ELT(val, 1, ssclme_sigma(x, REML));
    sigmasq = REAL(VECTOR_ELT(val, 1))[0];
    sigmasq = (sigmasq) * (sigmasq);
    for (i = 0; i < nf; i++) {
	double *mm = REAL(VECTOR_ELT(Omg, i));
	int j, k, nci = nc[i], ncip1 = nci+1;

	F77_CALL(dpotrf)("U", &nci, mm, &nci, &j);
	if (j)			/* shouldn't happen */
	    error("DPOTRF returned error code %d on Omega[%d]",
		  j, i + 1);
	F77_CALL(dpotri)("U", &nci, mm, &nci, &j);
	if (j)			/* shouldn't happen */
	    error("DTRTRI returned error code %d on Omega[%d]",
		  j, i + 1);
	for (j = 0; j < nci; j++) {
	    mm[j * ncip1] *= sigmasq;
	    for (k = 0; k < j; k++) {
		mm[j + k * nci] = (mm[k + j * nci] *= sigmasq);
	    }
	}
    }
    UNPROTECT(2);
    return val;
}