SCM

SCM Repository

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

Annotation of /pkg/src/sscCrosstab.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 149 - (view) (download) (as text)

1 : bates 10 #include "sscCrosstab.h"
2 :    
3 :     SEXP sscCrosstab(SEXP flist, SEXP upper)
4 :     {
5 : bates 146 int **fpt, *Ap, *Gp, *TTi, *Ti, *Tj, *dims, i,
6 : bates 10 ncol = 0,
7 :     nfac = length(flist),
8 :     nfc2 = (nfac * (nfac - 1))/2, /* nfac choose 2 */
9 :     nobs = length(VECTOR_ELT(flist, 0)),
10 : bates 146 ntrpl, nz, pos, up = asLogical(upper);
11 :     double *TTx, *Tx;
12 : bates 10 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("sscCrosstab")));
13 :    
14 :     if (!isNewList(flist) || nfac < 1)
15 :     error("flist must be a non-empty list");
16 :     SET_SLOT(val, Matrix_GpSym, allocVector(INTSXP, nfac + 1));
17 :     Gp = INTEGER(GET_SLOT(val, Matrix_GpSym));
18 :     fpt = (int **) R_alloc(nfac, sizeof(int *));
19 :     for (i = 0; i < nfac; i++) {
20 :     SEXP el = VECTOR_ELT(flist, i);
21 :     if (!inherits(el, "factor"))
22 :     error("flist must be a non-empty list of factors");
23 :     if (length(el) != nobs)
24 :     error("All elements of flist must have the same length");
25 :     Gp[i] = ncol;
26 :     ncol += length(getAttrib(el, R_LevelsSymbol));
27 :     fpt[i] = INTEGER(el);
28 :     }
29 :     Gp[nfac] = ncol;
30 :     SET_SLOT(val, Matrix_uploSym, ScalarString(mkChar(up ? "U" : "L")));
31 :     SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
32 :     dims = INTEGER(GET_SLOT(val, Matrix_DimSym));
33 :     dims[0] = dims[1] = ncol;
34 :     ntrpl = nfc2 * nobs + ncol;
35 :     Ti = Calloc(ntrpl, int); Tj = Calloc(ntrpl, int); TTi = Calloc(ntrpl, int);
36 :     Tx = Calloc(ntrpl, double); TTx = Calloc(ntrpl, double);
37 :     /* Generate the triplet form of the result */
38 :     for (i = 0; i < ncol; i++) {
39 :     Ti[i] = Tj[i] = i; /* The diagonals - these will store counts */
40 :     Tx[i] = 0.0;
41 :     }
42 :     pos = ncol;
43 :     for (i = 0; i < nobs ; i++) {
44 :     int j, jcol, k;
45 :     for (j = 0; j < nfac; j++) {
46 :     jcol = Gp[j] + fpt[j][i] - 1;
47 :     Tx[jcol] += 1.; /* increment diagonal count */
48 :     for (k = j + 1; k < nfac; k++) { /* off-diagonals */
49 :     int irow = Gp[k] + fpt[k][i] - 1;
50 :     if (up) {
51 :     Ti[pos] = jcol; Tj[pos] = irow;
52 :     } else {
53 :     Tj[pos] = jcol; Ti[pos] = irow;
54 :     }
55 :     Tx[pos] = 1.;
56 :     pos++;
57 :     }
58 :     }
59 :     }
60 :     SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, ncol + 1));
61 :     Ap = INTEGER(GET_SLOT(val, Matrix_pSym));
62 :     triplet_to_col(ncol, ncol, ntrpl, Ti, Tj, Tx, Ap, TTi, TTx);
63 :     nz = Ap[ncol]; /* non-zeros in Z'Z crosstab */
64 :     SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nz));
65 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nz));
66 :     Memcpy(INTEGER(GET_SLOT(val, Matrix_iSym)), TTi, nz);
67 :     Memcpy(REAL(GET_SLOT(val, Matrix_xSym)), TTx, nz);
68 :     Free(Ti); Free(Tj); Free(Tx);
69 :     Free(TTi); Free(TTx);
70 :    
71 :     UNPROTECT(1);
72 :     return val;
73 :     }
74 : bates 73
75 :     SEXP sscCrosstab_L_LI_sizes(SEXP ctab, SEXP permexp)
76 :     {
77 :     SEXP ans = PROTECT(allocVector(INTSXP, 4));
78 :     int *Ai = INTEGER(GET_SLOT(ctab, Matrix_iSym)),
79 :     *Ap = INTEGER(GET_SLOT(ctab, Matrix_pSym)),
80 :     *aa = INTEGER(ans),
81 :     *perm = INTEGER(permexp),
82 :     n = INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],
83 :     *Lp = Calloc(n + 1, int),
84 :     *Parent = Calloc(n, int),
85 :     *Lnz = Calloc(n, int),
86 :     *Flag = Calloc(n, int);
87 :    
88 :     ldl_symbolic(n, Ap, Ai, Lp, Parent, Lnz, Flag,
89 :     (int *) NULL, (int *) NULL); /* P & Pinv */
90 :     aa[0] = Lp[n];
91 :     ssclme_fill_LIp(n, Parent, Lp);
92 :     aa[1] = Lp[n];
93 :     ssc_symbolic_permute(n, 1, perm, Ap, Ai);
94 :     ldl_symbolic(n, Ap, Ai, Lp, Parent, Lnz, Flag,
95 :     (int *) NULL, (int *) NULL); /* P & Pinv */
96 :     aa[2] = Lp[n];
97 :     ssclme_fill_LIp(n, Parent, Lp);
98 :     aa[3] = Lp[n];
99 :     Free(Flag); Free(Lnz); Free(Parent); Free(Lp);
100 :     UNPROTECT(1);
101 :     return ans;
102 :     }
103 : bates 85
104 : bates 143 static
105 : bates 148 void col_metis_order(int j0, int j1, int i2,
106 :     const int Tp[], const int Ti[], int ans[])
107 : bates 143 {
108 : bates 148 int j, nz = 0; /* count off-diagonal pairs */
109 :     for (j = j0; j < j1; j++) { /* columns of interest */
110 :     int ii, nr = 0, p2 = Tp[j + 1];
111 :     for (ii = Tp[j]; ii < p2; ii++) {
112 :     int i = Ti[ii];
113 :     if (j1 <= i && i < i2) nr++; /* verify row index */
114 : bates 143 }
115 : bates 148 nz += (nr * (nr - 1))/2; /* add number of pairs of rows */
116 : bates 143 }
117 : bates 148 if (nz > 0) { /* Form an ssc Matrix */
118 :     int j, n = i2 - j1, /* number of rows */
119 :     nnz = n + nz, pos;
120 :     int *Ap = Calloc(n + 1, int),
121 :     *Ai = Calloc(nnz, int),
122 :     *Tj = Calloc(nnz, int),
123 : bates 149 *TTi = Calloc(nnz, int),
124 :     *perm = Calloc(n, int),
125 :     *iperm = Calloc(n, int);
126 : bates 143
127 : bates 148 for (j = 0; j < n; j++) { /* diagonals */
128 :     TTi[j] = Tj[j] = j;
129 :     Tx[j] = 1.;
130 : bates 85 }
131 : bates 148 pos = n;
132 :     for (j = j0; j < j1; j++) { /* create the pairs */
133 :     int ii, nr = 0, p2 = Tp[j + 1];
134 :     for (ii = Tp[j]; ii < p2; ii++) {
135 :     int r1 = Ti[ii], i1;
136 :     if (j1 <= r1 && r1 < i2) {
137 :     for (i1 = ii + 1; i1 < p2; i1++) {
138 :     int r2 = Ti[i1];
139 :     if (r2 < i2) {
140 :     TTi[pos] = r2 - j1;
141 :     Tj[pos] = r1 - j1;
142 :     Tx[pos] = 1.;
143 :     pos++;
144 :     }
145 :     }
146 : bates 93 }
147 : bates 85 }
148 :     }
149 : bates 149 triplet_to_col(n, n, nnz, TTi, Tj, (double *) 0, Ap, Ai, (double *) 0);
150 : bates 148 ssc_metis_order(n, Ap, Ai, perm, iperm);
151 :     for (j = j1; j < i2; j++) ans[j] = j1 + iperm[j - j1];
152 : bates 149 Free(TTi); Free(Tj); Free(Ai); Free(Ap);
153 : bates 148 Free(perm); Free(iperm);
154 : bates 85 }
155 :     }
156 :    
157 : bates 127 SEXP sscCrosstab_groupedPerm(SEXP ctab)
158 : bates 74 {
159 :     SEXP
160 :     GpSlot = GET_SLOT(ctab, Matrix_GpSym),
161 :     iSlot = GET_SLOT(ctab, Matrix_iSym),
162 :     pSlot = GET_SLOT(ctab, Matrix_pSym);
163 :     int *Ai = INTEGER(iSlot),
164 :     *Ap = INTEGER(pSlot),
165 :     *Gp = INTEGER(GpSlot),
166 : bates 85 i,
167 : bates 74 n = length(pSlot) - 1, /* number of columns */
168 :     nf = length(GpSlot) - 1, /* number of factors */
169 : bates 143 up;
170 : bates 74 SEXP ans = PROTECT(allocVector(INTSXP, n));
171 :    
172 : bates 143 up = toupper(*CHAR(STRING_ELT(GET_SLOT(ctab, Matrix_uploSym), 0))) != 'L';
173 :     if (nf > 1 && up) { /* transpose */
174 :     int nz = length(iSlot);
175 :     int *ai = Calloc(nz, int),
176 :     *ap = Calloc(n + 1, int);
177 :     double *ax = Calloc(nz, double);
178 : bates 74
179 : bates 143 csc_components_transpose(n, n, nz, Ap, Ai,
180 :     REAL(GET_SLOT(ctab, Matrix_xSym)),
181 :     ap, ai, ax);
182 :     Ap = ap;
183 :     Ai = ai;
184 :     Free(ax); /* don't need values, only positions */
185 :     }
186 : bates 85 for (i = 0; i < n; i++) {
187 : bates 86 INTEGER(ans)[i] = i; /* initialize permutation to identity */
188 : bates 85 }
189 : bates 146 for (i = 1; i < nf; i++) {
190 :     col_metis_order(Gp[i - 1], Gp[i], Gp[i+1], Ap, Ai, INTEGER(ans));
191 :     }
192 : bates 143 if (nf > 1 && up) {Free(Ap); Free(Ai);}
193 : bates 74 UNPROTECT(1);
194 :     return ans;
195 :     }

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