SCM

SCM Repository

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

Annotation of /pkg/src/lgCMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : bates 692 #include "lgCMatrix.h"
2 :    
3 :     SEXP lgCMatrix_validate(SEXP x)
4 :     {
5 :     SEXP pslot = GET_SLOT(x, Matrix_pSym),
6 :     islot = GET_SLOT(x, Matrix_iSym);
7 :     int j,
8 :     ncol = length(pslot) - 1,
9 :     *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
10 :     nrow,
11 :     *xp = INTEGER(pslot),
12 :     *xi = INTEGER(islot);
13 :    
14 :     nrow = dims[0];
15 :     if (length(pslot) <= 0)
16 :     return mkString(_("slot p must have length > 0"));
17 :     if (xp[0] != 0)
18 :     return mkString(_("first element of slot p must be zero"));
19 :     if (length(islot) != xp[ncol])
20 :     return mkString(_("last element of slot p must match length of slot i"));
21 :     for (j = 0; j < ncol; j++) {
22 :     if (xp[j] > xp[j+1])
23 :     return mkString(_("slot p must be non-decreasing"));
24 :     }
25 :     for (j = 0; j < length(islot); j++) {
26 :     if (xi[j] < 0 || xi[j] >= nrow)
27 :     return mkString(_("all row indices must be between 0 and nrow-1"));
28 :     }
29 :     if (csc_unsorted_columns(ncol, xp, xi)) {
30 :     csc_sort_columns(ncol, xp, xi, (double *) NULL);
31 :     }
32 :     return ScalarLogical(1);
33 :     }
34 :    
35 : maechler 866 /* very parallel to csc_matrix() in ./dgeMatrix.c */
36 :     SEXP lcsc_to_matrix(SEXP x)
37 :     {
38 :     SEXP ans, pslot = GET_SLOT(x, Matrix_pSym);
39 :     int j, ncol = length(pslot) - 1,
40 :     nrow = INTEGER(GET_SLOT(x, Matrix_DimSym))[0],
41 :     *xp = INTEGER(pslot),
42 :     *xi = INTEGER(GET_SLOT(x, Matrix_iSym));
43 :     int *ax;
44 :    
45 :     ax = LOGICAL(ans = PROTECT(allocMatrix(LGLSXP, nrow, ncol)));
46 :     for (j = 0; j < (nrow * ncol); j++) ax[j] = 0;
47 :     for (j = 0; j < ncol; j++) {
48 :     int ind;
49 :     for (ind = xp[j]; ind < xp[j+1]; ind++)
50 :     ax[j * nrow + xi[ind]] = 1;
51 :     }
52 :     UNPROTECT(1);
53 :     return ans;
54 :     }
55 :    
56 :     /**
57 : bates 692 * C := op(A) %*% op(B) + beta ^ C for logical sparse column-oriented matrices
58 : maechler 866 *
59 : bates 692 * @param tra nonzero if A is to be transposed
60 :     * @param trb nonzero if B is to be transposed
61 :     * @param m number of rows in C
62 :     * @param n number of columns in C
63 :     * @param k number of columns in A if tra == 0, otherwise number of
64 :     * rows in A
65 :     * @param ai vector of row indices of TRUE elements in A
66 :     * @param ap column pointers for A
67 :     * @param bi vector of row indices of TRUE elements in B
68 :     * @param bp column pointers for B
69 :     * @param beta if non-zero existing TRUE elements in C are retained
70 : bates 718 * @param ciP SEXP whose INTEGER part is the column indices of TRUE
71 :     * elements in C (not used if beta == 0).
72 : bates 692 * @param cp column pointers for C
73 : bates 718 *
74 :     * @return SEXP whose INTEGER part is the column indices of TRUE
75 :     * elements in the product. Note that the contents of cp may be modified.
76 : bates 692 */
77 : bates 718 SEXP Matrix_lgClgCmm(int tra, int trb, int m, int n, int k,
78 : bates 692 const int ai[], const int ap[],
79 :     const int bi[], const int bp[],
80 : bates 718 int beta, SEXP CIP, int cp[])
81 : bates 692 {
82 : bates 718 int cnnz = cp[n], extra = 0;
83 :     int *ci, i, j, prot = 0; /* prot is the number of PROTECTs to UNPROTECT */
84 : maechler 866
85 : bates 692 if (beta) {
86 : bates 718 ci = INTEGER(CIP);
87 :     } else { /* blank the C matrix */
88 :     for (j = 0; j <= n; j++) cp[j] = 0;
89 :     cnnz = 0;
90 :     ci = (int *) NULL;
91 : maechler 866 }
92 : bates 723
93 :     if (tra) { /* replace ai and ap by els for transpose */
94 :     int nz = ap[m];
95 :     int *Ai = Calloc(nz, int),
96 : bates 732 *aj = expand_cmprPt(m, ap, Calloc(nz, int)),
97 : bates 723 *Ap = Calloc(k + 1, int);
98 : maechler 866
99 : bates 732 triplet_to_col(m, k, nz, aj, ai, (double *) NULL,
100 :     Ap, Ai, (double *) NULL);
101 :     Free(aj);
102 : bates 723 ai = Ai; ap = Ap;
103 :     }
104 :    
105 :     if (trb) { /* replace bi and bp by els for transpose */
106 :     int nz = bp[k];
107 :     int *Bi = Calloc(nz, int),
108 : bates 732 *bj = expand_cmprPt(k, bp, Calloc(nz, int)),
109 : bates 723 *Bp = Calloc(n + 1, int);
110 : maechler 866
111 : bates 732 triplet_to_col(k, n, nz, bj, bi, (double *) NULL,
112 :     Bp, Bi, (double *) NULL);
113 :     Free(bj);
114 : bates 723 bi = Bi; bp = Bp;
115 :     }
116 :    
117 :     for (j = 0; j < n; j++) { /* col index for B and C */
118 :     int ii, ii2 = bp[j + 1];
119 :     for (ii = bp[j]; ii < ii2; ii++) { /* index into bi */
120 :     int jj = bi[ii]; /* row index of B; col index of A */
121 :     int i, i2 = ap[jj + 1]; /* index into ai */
122 :     for (i = ap[jj]; i < i2; i++)
123 :     if (check_csc_index(cp, ci, ai[i], j, -1) < 0) extra++;
124 : bates 718 }
125 : bates 723 }
126 : bates 718
127 : bates 723 if (extra) {
128 :     int ntot = cnnz + extra;
129 : bates 732 int *Cp = Calloc(n + 1, int),
130 :     *Ti = Calloc(ntot, int),
131 : bates 723 *rwInd = Calloc(m, int), /* indicator of TRUE in column j */
132 :     pos = 0;
133 : maechler 866
134 : bates 732 Cp[0] = 0;
135 : bates 723 for (j = 0; j < n; j++) {
136 :     int ii, ii2 = bp[j + 1];
137 : maechler 866
138 : bates 732 AZERO(rwInd, m); /* initialize column j of C */
139 :     for (i = cp[j]; i < cp[j+1]; i++) rwInd[ci[i]] = 1;
140 :    
141 :     Cp[j + 1] = Cp[j];
142 : bates 723 for (ii = bp[j]; ii < ii2; ii++) { /* index into bi */
143 :     int jj = bi[ii]; /* row index of B; col index of A */
144 :     int i, i2 = ap[jj + 1]; /* index into ai */
145 :     for (i = ap[jj]; i < i2; i++) rwInd[ai[i]] = 1;
146 : bates 692 }
147 : bates 723 for (i = 0; i < m; i++)
148 : bates 732 if (rwInd[i]) {Cp[j + 1]++; Ti[pos++] = i;}
149 : bates 692 }
150 : bates 732 PROTECT(CIP = allocVector(INTSXP, Cp[n])); prot++;
151 :     Memcpy(INTEGER(CIP), Ti, Cp[n]);
152 :     Memcpy(cp, Cp, n + 1);
153 :     Free(Cp); Free(Ti); Free(rwInd);
154 : bates 692 }
155 : bates 723
156 :     if (tra) {Free(ai); Free(ap);}
157 :     if (trb) {Free(bi); Free(bp);}
158 : bates 718 UNPROTECT(prot);
159 :     return CIP;
160 : bates 692 }
161 : maechler 866
162 : bates 718 SEXP lgCMatrix_lgCMatrix_mm(SEXP a, SEXP b)
163 :     {
164 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lgCMatrix")));
165 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
166 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
167 :     *cdims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
168 :     int k = adims[1], m = adims[0], n = bdims[1];
169 :     int *cp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
170 : maechler 866
171 : bates 718 if (bdims[0] != k)
172 :     error(_("Matrices are not conformable for multiplication"));
173 :     cdims[0] = m; cdims[1] = n;
174 :     SET_SLOT(ans, Matrix_iSym,
175 :     Matrix_lgClgCmm(0, 0, m, n, k,
176 :     INTEGER(GET_SLOT(a, Matrix_iSym)),
177 :     INTEGER(GET_SLOT(a, Matrix_pSym)),
178 :     INTEGER(GET_SLOT(b, Matrix_iSym)),
179 :     INTEGER(GET_SLOT(b, Matrix_pSym)),
180 :     0, (SEXP) NULL, cp));
181 :     UNPROTECT(1);
182 :     return ans;
183 :     }
184 : bates 692
185 : bates 718 SEXP lgCMatrix_trans(SEXP x)
186 : bates 692 {
187 : bates 718 SEXP xi = GET_SLOT(x, Matrix_iSym);
188 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lgCMatrix")));
189 :     int *adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)),
190 :     *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
191 :     nz = length(xi);
192 :     int *xj = Calloc(nz, int);
193 : bates 726 SEXP adn = ALLOC_SLOT(ans, Matrix_DimNamesSym, VECSXP, 2),
194 :     xdn = GET_SLOT(x, Matrix_DimNamesSym);
195 : maechler 866
196 : bates 726 adims[1] = xdims[0]; adims[0] = xdims[1];
197 :     SET_VECTOR_ELT(adn, 0, VECTOR_ELT(xdn, 1));
198 :     SET_VECTOR_ELT(adn, 1, VECTOR_ELT(xdn, 0));
199 : maechler 866 triplet_to_col(adims[0], adims[1], nz,
200 : bates 718 expand_cmprPt(xdims[1], INTEGER(GET_SLOT(x, Matrix_pSym)), xj),
201 : bates 726 INTEGER(xi), (double *) NULL,
202 : bates 718 INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, adims[1] + 1)),
203 :     INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)),
204 :     (double *) NULL);
205 :     Free(xj);
206 :     UNPROTECT(1);
207 :     return ans;
208 : bates 692 }
209 : bates 723
210 : maechler 866 /**
211 : bates 723 * Replace C by AA' + beta*C or A'A + beta*C
212 : maechler 866 *
213 : bates 723 * @param up Indicator of upper/lower triangle in the symmetric sparse matrix
214 :     * @param tra Transpose, in the sense of dsyrk. That is, tra TRUE indicates A'A
215 :     * @param n size of the product matrix
216 :     * @param k number of columns in A if tra is FALSE, otherwise the number of rows
217 :     * @param ai row indices for A
218 :     * @param ap column pointers for A
219 :     * @param beta TRUE if existing elements in C are to be preserved
220 :     * @param CIP SEXP whose INTEGER part is the row indices of C (not used if beta is FALSE)
221 :     * @param cp column pointers for C
222 : maechler 866 *
223 : bates 723 * @return SEXP whose INTEGER part is the updated row indices of C
224 :     */
225 :     SEXP Matrix_lgCsyrk(int up, int tra, int n, int k, const int ai[], const int ap[],
226 :     int beta, SEXP CIP, int cp[])
227 :     {
228 :     int extra = 0, i, ii, j, prot = 0;
229 :     int *ci, cnnz = cp[n];
230 :    
231 :     if (beta) {
232 :     ci = INTEGER(CIP);
233 :     } else { /* blank the C matrix */
234 :     for (j = 0; j <= n; j++) cp[j] = 0;
235 :     cnnz = 0;
236 :     ci = (int *) NULL;
237 : maechler 866 }
238 :    
239 : bates 723 if (tra) { /* replace ai and ap by els for transpose */
240 :     int nz = ap[n];
241 :     int *Ai = Calloc(nz, int),
242 : bates 732 *aj = expand_cmprPt(n, ap, Calloc(nz, int)),
243 : bates 723 *Ap = Calloc(k + 1, int);
244 : maechler 866
245 : bates 732 triplet_to_col(n, k, nz, aj, ai, (double *) NULL,
246 :     Ap, Ai, (double *) NULL);
247 :     Free(aj);
248 : bates 723 ai = Ai; ap = Ap;
249 :     }
250 : bates 732
251 : bates 723 for (j = 0; j < k; j++) {
252 :     int i2 = ap[j + 1];
253 :     for (i = ap[j]; i < i2; i++) {
254 :     int r1 = ai[i];
255 : bates 732 if (r1 < 0 || r1 >= n)
256 :     error(_("row %d not in row range [0,%d]"), r1, n - 1);
257 : bates 723 for (ii = i; ii < i2; ii++) {
258 :     int r2 = ai[ii];
259 : bates 732 if (r2 < 0 || r2 >= n)
260 :     error(_("row %d not in row range [0,%d]"), r2, n - 1);
261 : bates 723 if (check_csc_index(cp, ci, up?r1:r2, up?r2:r1, -1) < 0)
262 :     extra++;
263 :     }
264 :     }
265 :     }
266 :    
267 :     if (extra) {
268 :     int ntot = cnnz + extra;
269 :     int *Ti = Memcpy(Calloc(ntot, int), ci, cnnz),
270 :     *Tj = expand_cmprPt(n, cp, Calloc(ntot, int)),
271 :     *Ci = Calloc(ntot, int),
272 :     pos = cnnz;
273 : maechler 866
274 : bates 723 for (j = 0; j < k; j++) {
275 :     int i2 = ap[j + 1];
276 :     for (i = ap[j]; i < i2; i++) {
277 :     int r1 = ai[i];
278 :     for (ii = i; ii < i2; ii++) {
279 :     int r2 = ai[ii];
280 :     int row = up ? r1 : r2, col = up ? r2 : r1;
281 :     if (r2 < r1) error("[j,i,ii,r1,r2] = [%d,%d,%d,%d,%d]",
282 :     j,i,ii,r1,r2);
283 :     if (check_csc_index(cp, ci, row, col, -1) < 0) {
284 :     Ti[pos] = row;
285 :     Tj[pos] = col;
286 :     pos++;
287 :     }
288 :     }
289 :     }
290 :     }
291 : maechler 866
292 : bates 723 triplet_to_col(n, n, pos, Ti, Tj, (double *) NULL,
293 :     cp, Ci, (double *) NULL);
294 :     PROTECT(CIP = allocVector(INTSXP, cp[n])); prot++;
295 :     Memcpy(INTEGER(CIP), Ci, cp[n]);
296 :     Free(Ti); Free(Tj); Free(Ci);
297 :     }
298 : maechler 866
299 : bates 723 if (tra) {Free(ai); Free(ap);}
300 :     UNPROTECT(prot);
301 :     return CIP;
302 :     }
303 :    
304 : maechler 866 /**
305 : bates 723 * Create the cross-product or transpose cross-product of a logical
306 :     * sparse matrix in column-oriented compressed storage mode.
307 : maechler 866 *
308 : bates 723 * @param x Pointer to a lgCMatrix
309 :     * @param trans logical indicator of transpose, in the sense of dsyrk.
310 : bates 732 * That is, trans == TRUE is used for crossprod.
311 : maechler 866 * @param C
312 :     *
313 : bates 723 * @return An lsCMatrix of the form if(trans) X'X else XX'
314 :     */
315 : bates 732 SEXP lgCMatrix_crossprod(SEXP x, SEXP trans, SEXP C)
316 : bates 723 {
317 :     int tra = asLogical(trans);
318 : bates 732 int *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
319 : bates 723 int k = xdims[tra ? 0 : 1], n = xdims[tra ? 1 : 0];
320 : maechler 866
321 : bates 732 if (C == R_NilValue) {
322 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lsCMatrix")));
323 :    
324 :     adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
325 :     adims[0] = adims[1] = n;
326 :     SET_SLOT(ans, Matrix_uploSym, mkString("U"));
327 :     SET_SLOT(ans, Matrix_iSym,
328 :     Matrix_lgCsyrk(1, tra, n, k,
329 :     INTEGER(GET_SLOT(x, Matrix_iSym)),
330 :     INTEGER(GET_SLOT(x, Matrix_pSym)),
331 :     0, R_NilValue,
332 :     INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1))));
333 :     UNPROTECT(1);
334 :     return ans;
335 :     }
336 :     adims = INTEGER(GET_SLOT(C, Matrix_DimSym));
337 :     if (adims[0] != n || adims[1] != n)
338 :     error(_("Dimensions of x and y are not compatible for crossprod"));
339 :     SET_SLOT(C, Matrix_iSym,
340 :     Matrix_lgCsyrk(CHAR(asChar(GET_SLOT(C, Matrix_uploSym)))[0] == 'U',
341 :     tra, n, k,
342 : bates 723 INTEGER(GET_SLOT(x, Matrix_iSym)),
343 :     INTEGER(GET_SLOT(x, Matrix_pSym)),
344 : bates 732 1, GET_SLOT(C, Matrix_iSym),
345 :     INTEGER(GET_SLOT(C, Matrix_pSym))));
346 :     return C;
347 :     }
348 :    
349 : maechler 866 /**
350 : bates 732 * Special-purpose function that returns a permutation of the columns
351 :     * of a lgTMatrix for which nrow(x) > ncol(x). The ordering puts
352 :     * columns with fewer entries on the left. Once a column has been
353 :     * moved to the left the rows in where that column is TRUE are removed
354 :     * from the counts.
355 : maechler 866 *
356 : bates 732 * @param x Pointer to an lgTMatrix object
357 : maechler 866 *
358 : bates 732 * @return 0-based permutation vector for the columns of x
359 :     */
360 :     SEXP lgCMatrix_picky_column(SEXP x)
361 :     {
362 :     int *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
363 :     int *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
364 :     *xp = INTEGER(GET_SLOT(x, Matrix_pSym)),
365 :     m = xdims[0], n = xdims[1];
366 :     SEXP ans = PROTECT(allocVector(INTSXP, n));
367 :     int *actr = Calloc(m, int),
368 :     *actc = Calloc(n, int),
369 : bates 762 cj, i, j, mincount, minloc = -1, pos;
370 : bates 732
371 :     for (i = 0; i < m; i++) actr[i] = 1;
372 :     mincount = m + 1;
373 :     for (j = 0; j < n; j++) {
374 :     cj = xp[j + 1] - xp[j];
375 :     actc[j] = 1;
376 :     if (cj < mincount) {
377 :     mincount = cj;
378 :     minloc = j;
379 :     }
380 :     }
381 : maechler 866
382 : bates 732 pos = 0;
383 :     while (pos < n) {
384 :     INTEGER(ans)[pos++] = minloc;
385 :     actc[minloc] = 0;
386 :     for (i = xp[minloc]; i < xp[minloc + 1]; i++) actr[xi[i]] = 0;
387 :     mincount = m + 1;
388 :     for (j = 0; j < n; j++) {
389 :     if (actc[j]) {
390 :     cj = 0;
391 :     for (i = xp[j]; i < xp[j + 1]; i++) {
392 :     if (actr[xi[i]]) cj++;
393 :     if (cj < mincount) {
394 :     mincount = cj;
395 :     minloc = j;
396 :     }
397 :     }
398 :     }
399 :     }
400 :     }
401 :    
402 :     Free(actr); Free(actc);
403 : bates 723 UNPROTECT(1);
404 :     return ans;
405 :     }

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