# SCM Repository

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

# Diff of /pkg/src/ssclme.c

revision 125, Sun Apr 25 11:59:37 2004 UTC revision 143, Fri Apr 30 20:33:01 2004 UTC
# Line 1  Line 1
1  #include "ssclme.h"  #include "ssclme.h"
2
/**
* 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;
}

3  static  static
4  void ssclme_copy_ctab(int nf, const int nc[], SEXP ctab, SEXP ssc)  void ssclme_copy_ctab(int nf, const int nc[], SEXP ctab, SEXP ssc)
5  {  {
# Line 230  Line 141
141      ssc = VECTOR_ELT(val, 0);      ssc = VECTOR_ELT(val, 0);
142                                  /* Pairwise cross-tabulation */                                  /* Pairwise cross-tabulation */
143      ctab = PROTECT(sscCrosstab(facs, ScalarLogical(1)));      ctab = PROTECT(sscCrosstab(facs, ScalarLogical(1)));
144      SET_VECTOR_ELT(val, 1, ctab_permute(ctab));      SET_VECTOR_ELT(val, 1, sscCrosstab_groupedPerm(ctab));
145      if (length(VECTOR_ELT(val, 1)) > 0) {/* Fill-reducing permutation */      if (length(VECTOR_ELT(val, 1)) > 0) {/* Fill-reducing permutation */
146          ssc_symbolic_permute(INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],          ssc_symbolic_permute(INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],
147                               1, INTEGER(VECTOR_ELT(val, 1)),                               1, INTEGER(VECTOR_ELT(val, 1)),
# Line 369  Line 280
280          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
281          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
282          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
283            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
284          i, j, k,          i, j, k,
285          ione = 1,          ione = 1,
286          nf = length(mmats) - 1,          nf = length(mmats) - 1,
# Line 471  Line 383
383      }      }
384      Free(Z);      Free(Z);
385      ssclme_transfer_dimnames(x, facs, mmats);      ssclme_transfer_dimnames(x, facs, mmats);
386        status[0] = status[1] = 0;
387      return R_NilValue;      return R_NilValue;
388  }  }
389
# Line 1263  Line 1176
1176          alpha,          alpha,
1177          one = 1.;          one = 1.;
1178      SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));      SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
1179
1180      ssclme_factor(x);      ssclme_factor(x);
1181      if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0])) {      if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0])) {

Legend:
 Removed from v.125 changed lines Added in v.143