SCM

SCM Repository

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

Diff of /pkg/src/sscCrosstab.c

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

revision 183, Sun May 30 07:36:30 2004 UTC revision 184, Sun May 30 22:38:37 2004 UTC
# Line 170  Line 170 
170   *   *
171   * @param ctab pointer to a sscCrosstab object   * @param ctab pointer to a sscCrosstab object
172   *   *
173   * @return a pointer to an tscMatrix giving the projection of the 2,1 component   * @return a pointer to an sscMatrix giving the projection of the 2,1 component
174   */   */
175  SEXP sscCrosstab_project(SEXP ctab)  SEXP sscCrosstab_project(SEXP ctab)
176  {  {
# Line 261  Line 261 
261      UNPROTECT(1);      UNPROTECT(1);
262      return ans;      return ans;
263  }  }
264    
265    /**
266     * Project the first group of columns in an sscCrosstab object onto the
267     * remaining columns.
268     *
269     * @param ctab pointer to a sscCrosstab object
270     *
271     * @return a pointer to an sscMatrix with the projection
272     */
273    SEXP sscCrosstab_project2(SEXP ctab)
274    {
275        SEXP
276            GpSlot = GET_SLOT(ctab, Matrix_GpSym),
277            iSlot = GET_SLOT(ctab, Matrix_iSym),
278            pSlot = GET_SLOT(ctab, Matrix_pSym);
279        int *Ai = INTEGER(iSlot),
280            *Ap = INTEGER(pSlot),
281            *Gp = INTEGER(GpSlot),
282            i, i2, j, j1, k, k2,
283            nf = length(GpSlot) - 1, /* number of factors */
284            nz, pos, up, *AAp, *Ti, *Cp, *ind;
285        double *Ax = REAL(GET_SLOT(ctab, Matrix_xSym));
286        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sscMatrix")));
287    
288    
289        if (nf < 2) error("sscCrosstab_project2 requires more than one group");
290        up = toupper(*CHAR(STRING_ELT(GET_SLOT(ctab, Matrix_uploSym), 0))) != 'L';
291        if (up) {                   /* tranpose */
292            int n = length(pSlot) - 1;      /* number of columns */
293            int nnz = length(iSlot);
294            int *ai = Calloc(nnz, int);
295            int *ap = Calloc(n + 1, int);
296            double *ax = Calloc(nnz, double);
297    
298            csc_components_transpose(n, n, nnz, Ap, Ai, Ax, ap, ai, ax);
299            Ap = ap; Ai = ai; Ax = ax;
300        }
301    
302        j1 = Gp[1]; i2 = Gp[nf];    /* group boundaries */
303        Cp = Calloc(j1, int);       /* first pos in col with row > i */
304        for (j = 0; j < j1; j++) Cp[j] = Ap[j] + 1; /* Ap[j] is the diagonal */
305    
306        nz = Ap[i2] - Ap[j1];       /* upper bound on nonzeros in result */
307        for (j = 0; j < j1; j++) {
308            int nr = Ap[j + 1] - Ap[j] - 1;
309            nz += (nr * (nr - 1))/2; /* add number of pairs of rows below diag */
310        }
311    
312        Ti = Calloc(nz, int);       /* temporary row indices */
313        SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, i2 - j1 + 1));
314        AAp = INTEGER(GET_SLOT(ans, Matrix_pSym)); /* column pointers */
315    
316        AAp[0] = 0; pos = 0;
317        ind = Calloc(i2 - j1, int); /* indicator of rows in same column */
318        for (i = j1; i < i2; i++) {
319            for (k = j1; k < i2; k++) ind[k - j1] = 0;
320            for (j = 0; j < j1; j++) {
321                if (Ai[Cp[j]] == i) { /* go down the column */
322                    k2 = Ap[j+1];
323                    for(k = Cp[j] + 1; k < k2; k++) ind[Ai[k] - j1] = 1;
324                    Cp[j]++;
325                }
326            }
327    
328            Ti[pos++] = i - j1;     /* diagonal element */
329            for (k = i+1; k < i2; k++) { /* projected pairs */
330                int ii = k - j1;
331                if (ind[ii]) Ti[pos++] = ii;
332            }
333            k2 = Ap[i+1];
334            for (k = Ap[i] + 1; k < k2; k++) { /* previous off-diagonals */
335                Ti[pos++] = Ai[k] - j1;
336            }
337            AAp[i - j1 + 1] = pos;
338        }
339        nz = AAp[i2 - j1];
340        SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nz));
341        Memcpy(INTEGER(GET_SLOT(ans, Matrix_iSym)), Ti, nz);
342        SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nz));
343        Ax = REAL(GET_SLOT(ans, Matrix_xSym));
344        for (j = 0; j < nz; j++) Ax[j] = 1.;
345        SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));
346        SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
347        AAp = INTEGER(GET_SLOT(ans, Matrix_DimSym));
348        AAp[0] = AAp[1] = i2 - j1;
349        Free(Ti); Free(Cp); Free(ind);
350        if (up) {Free(Ap); Free(Ai); free(Ax);}
351        UNPROTECT(1);
352        return ans;
353    }

Legend:
Removed from v.183  
changed lines
  Added in v.184

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