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 723 - (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 :     /**
36 :     * C := op(A) %*% op(B) + beta ^ C for logical sparse column-oriented matrices
37 :     *
38 :     * @param tra nonzero if A is to be transposed
39 :     * @param trb nonzero if B is to be transposed
40 :     * @param m number of rows in C
41 :     * @param n number of columns in C
42 :     * @param k number of columns in A if tra == 0, otherwise number of
43 :     * rows in A
44 :     * @param ai vector of row indices of TRUE elements in A
45 :     * @param ap column pointers for A
46 :     * @param bi vector of row indices of TRUE elements in B
47 :     * @param bp column pointers for B
48 :     * @param beta if non-zero existing TRUE elements in C are retained
49 : bates 718 * @param ciP SEXP whose INTEGER part is the column indices of TRUE
50 :     * elements in C (not used if beta == 0).
51 : bates 692 * @param cp column pointers for C
52 : bates 718 *
53 :     * @return SEXP whose INTEGER part is the column indices of TRUE
54 :     * elements in the product. Note that the contents of cp may be modified.
55 : bates 692 */
56 : bates 718 SEXP Matrix_lgClgCmm(int tra, int trb, int m, int n, int k,
57 : bates 692 const int ai[], const int ap[],
58 :     const int bi[], const int bp[],
59 : bates 718 int beta, SEXP CIP, int cp[])
60 : bates 692 {
61 : bates 718 int cnnz = cp[n], extra = 0;
62 :     int *ci, i, j, prot = 0; /* prot is the number of PROTECTs to UNPROTECT */
63 : bates 692
64 :     if (beta) {
65 : bates 718 ci = INTEGER(CIP);
66 :     } else { /* blank the C matrix */
67 :     for (j = 0; j <= n; j++) cp[j] = 0;
68 :     cnnz = 0;
69 :     ci = (int *) NULL;
70 :     }
71 : bates 723
72 :     if (tra) { /* replace ai and ap by els for transpose */
73 :     int nz = ap[m];
74 :     int *Ai = Calloc(nz, int),
75 :     *Ti = Calloc(nz, int),
76 :     *Ap = Calloc(k + 1, int);
77 :    
78 :     triplet_to_col(k, m, nz, expand_cmprPt(m, ap, Ti), ai,
79 :     (double *) NULL, Ap, Ai, (double *) NULL);
80 :     Free(Ti);
81 :     ai = Ai; ap = Ap;
82 :     }
83 :    
84 :     if (trb) { /* replace bi and bp by els for transpose */
85 :     int nz = bp[k];
86 :     int *Bi = Calloc(nz, int),
87 :     *Ti = Calloc(nz, int),
88 :     *Bp = Calloc(n + 1, int);
89 :    
90 :     triplet_to_col(n, k, nz, expand_cmprPt(k, bp, Ti), bi,
91 :     (double *) NULL, Bp, Bi, (double *) NULL);
92 :     Free(Ti);
93 :     bi = Bi; bp = Bp;
94 :     }
95 :    
96 :     for (j = 0; j < n; j++) { /* col index for B and C */
97 :     int ii, ii2 = bp[j + 1];
98 :     for (ii = bp[j]; ii < ii2; ii++) { /* index into bi */
99 :     int jj = bi[ii]; /* row index of B; col index of A */
100 :     int i, i2 = ap[jj + 1]; /* index into ai */
101 :     for (i = ap[jj]; i < i2; i++)
102 :     if (check_csc_index(cp, ci, ai[i], j, -1) < 0) extra++;
103 : bates 718 }
104 : bates 723 }
105 : bates 718
106 : bates 723 if (extra) {
107 :     int ntot = cnnz + extra;
108 :     int *Ti = Calloc(ntot, int),
109 :     *rwInd = Calloc(m, int), /* indicator of TRUE in column j */
110 :     pos = 0;
111 :    
112 :     cp[0] = 0;
113 :     for (j = 0; j < n; j++) {
114 :     int ii, ii2 = bp[j + 1];
115 :     cp[j + 1] = cp[j];
116 :     AZERO(rwInd, m); /* initialize column j of C to FALSE */
117 :     for (ii = bp[j]; ii < ii2; ii++) { /* index into bi */
118 :     int jj = bi[ii]; /* row index of B; col index of A */
119 :     int i, i2 = ap[jj + 1]; /* index into ai */
120 :     for (i = ap[jj]; i < i2; i++) rwInd[ai[i]] = 1;
121 : bates 692 }
122 : bates 723 for (i = 0; i < m; i++)
123 :     if (rwInd[i]) {cp[j + 1]++; Ti[pos++] = i;}
124 : bates 692 }
125 : bates 723 PROTECT(CIP = allocVector(INTSXP, cp[n])); prot++;
126 :     Memcpy(INTEGER(CIP), Ti, cp[n]);
127 :     Free(Ti); Free(rwInd);
128 : bates 692 }
129 : bates 723
130 :     if (tra) {Free(ai); Free(ap);}
131 :     if (trb) {Free(bi); Free(bp);}
132 : bates 718 UNPROTECT(prot);
133 :     return CIP;
134 : bates 692 }
135 : bates 718
136 :     SEXP lgCMatrix_lgCMatrix_mm(SEXP a, SEXP b)
137 :     {
138 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lgCMatrix")));
139 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
140 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
141 :     *cdims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
142 :     int k = adims[1], m = adims[0], n = bdims[1];
143 :     int *cp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
144 :    
145 :     if (bdims[0] != k)
146 :     error(_("Matrices are not conformable for multiplication"));
147 :     cdims[0] = m; cdims[1] = n;
148 :     SET_SLOT(ans, Matrix_iSym,
149 :     Matrix_lgClgCmm(0, 0, m, n, k,
150 :     INTEGER(GET_SLOT(a, Matrix_iSym)),
151 :     INTEGER(GET_SLOT(a, Matrix_pSym)),
152 :     INTEGER(GET_SLOT(b, Matrix_iSym)),
153 :     INTEGER(GET_SLOT(b, Matrix_pSym)),
154 :     0, (SEXP) NULL, cp));
155 :     UNPROTECT(1);
156 :     return ans;
157 :     }
158 : bates 692
159 : bates 718 SEXP lgCMatrix_trans(SEXP x)
160 : bates 692 {
161 : bates 718 SEXP xi = GET_SLOT(x, Matrix_iSym);
162 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lgCMatrix")));
163 :     int *adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)),
164 :     *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
165 :     nz = length(xi);
166 :     int *xj = Calloc(nz, int);
167 :    
168 :     adims[1] = xdims[0];
169 :     adims[0] = xdims[1];
170 :     triplet_to_col(adims[0], adims[1], nz,
171 :     expand_cmprPt(xdims[1], INTEGER(GET_SLOT(x, Matrix_pSym)), xj),
172 :     INTEGER(GET_SLOT(x, Matrix_iSym)), (double *) NULL,
173 :     INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, adims[1] + 1)),
174 :     INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)),
175 :     (double *) NULL);
176 :     Free(xj);
177 :     UNPROTECT(1);
178 :     return ans;
179 : bates 692 }
180 : bates 723
181 :     /**
182 :     * Replace C by AA' + beta*C or A'A + beta*C
183 :     *
184 :     * @param up Indicator of upper/lower triangle in the symmetric sparse matrix
185 :     * @param tra Transpose, in the sense of dsyrk. That is, tra TRUE indicates A'A
186 :     * @param n size of the product matrix
187 :     * @param k number of columns in A if tra is FALSE, otherwise the number of rows
188 :     * @param ai row indices for A
189 :     * @param ap column pointers for A
190 :     * @param beta TRUE if existing elements in C are to be preserved
191 :     * @param CIP SEXP whose INTEGER part is the row indices of C (not used if beta is FALSE)
192 :     * @param cp column pointers for C
193 :     *
194 :     * @return SEXP whose INTEGER part is the updated row indices of C
195 :     */
196 :     SEXP Matrix_lgCsyrk(int up, int tra, int n, int k, const int ai[], const int ap[],
197 :     int beta, SEXP CIP, int cp[])
198 :     {
199 :     int extra = 0, i, ii, j, prot = 0;
200 :     int *ci, cnnz = cp[n];
201 :    
202 :     if (beta) {
203 :     ci = INTEGER(CIP);
204 :     } else { /* blank the C matrix */
205 :     for (j = 0; j <= n; j++) cp[j] = 0;
206 :     cnnz = 0;
207 :     ci = (int *) NULL;
208 :     }
209 :    
210 :     if (tra) { /* replace ai and ap by els for transpose */
211 :     int nz = ap[n];
212 :     int *Ai = Calloc(nz, int),
213 :     *Ti = Calloc(nz, int),
214 :     *Ap = Calloc(k + 1, int);
215 :    
216 :     triplet_to_col(k, n, nz, expand_cmprPt(k, ap, Ti), ai,
217 :     (double *) NULL, Ap, Ai, (double *) NULL);
218 :     Free(Ti);
219 :     ai = Ai; ap = Ap;
220 :     }
221 :     for (j = 0; j < k; j++) {
222 :     int i2 = ap[j + 1];
223 :     for (i = ap[j]; i < i2; i++) {
224 :     int r1 = ai[i];
225 :     for (ii = i; ii < i2; ii++) {
226 :     int r2 = ai[ii];
227 :     if (check_csc_index(cp, ci, up?r1:r2, up?r2:r1, -1) < 0)
228 :     extra++;
229 :     }
230 :     }
231 :     }
232 :    
233 :     if (extra) {
234 :     int ntot = cnnz + extra;
235 :     int *Ti = Memcpy(Calloc(ntot, int), ci, cnnz),
236 :     *Tj = expand_cmprPt(n, cp, Calloc(ntot, int)),
237 :     *Ci = Calloc(ntot, int),
238 :     pos = cnnz;
239 :    
240 :     for (j = 0; j < k; j++) {
241 :     int i2 = ap[j + 1];
242 :     for (i = ap[j]; i < i2; i++) {
243 :     int r1 = ai[i];
244 :     for (ii = i; ii < i2; ii++) {
245 :     int r2 = ai[ii];
246 :     int row = up ? r1 : r2, col = up ? r2 : r1;
247 :     if (r2 < r1) error("[j,i,ii,r1,r2] = [%d,%d,%d,%d,%d]",
248 :     j,i,ii,r1,r2);
249 :     if (check_csc_index(cp, ci, row, col, -1) < 0) {
250 :     Ti[pos] = row;
251 :     Tj[pos] = col;
252 :     pos++;
253 :     }
254 :     }
255 :     }
256 :     }
257 :    
258 :     triplet_to_col(n, n, pos, Ti, Tj, (double *) NULL,
259 :     cp, Ci, (double *) NULL);
260 :     PROTECT(CIP = allocVector(INTSXP, cp[n])); prot++;
261 :     Memcpy(INTEGER(CIP), Ci, cp[n]);
262 :     Free(Ti); Free(Tj); Free(Ci);
263 :     }
264 :    
265 :     if (tra) {Free(ai); Free(ap);}
266 :     UNPROTECT(prot);
267 :     return CIP;
268 :     }
269 :    
270 :     /**
271 :     * Create the cross-product or transpose cross-product of a logical
272 :     * sparse matrix in column-oriented compressed storage mode.
273 :     *
274 :     * @param x Pointer to a lgCMatrix
275 :     * @param trans logical indicator of transpose, in the sense of dsyrk.
276 :     * That is, trans TRUE is used for crossprod.
277 :     *
278 :     * @return An lsCMatrix of the form if(trans) X'X else XX'
279 :     */
280 :     SEXP lgCMatrix_crossprod(SEXP x, SEXP trans)
281 :     {
282 :     int tra = asLogical(trans);
283 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lsCMatrix")));
284 :     int *adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)),
285 :     *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
286 :     int k = xdims[tra ? 0 : 1], n = xdims[tra ? 1 : 0];
287 :    
288 :     adims[0] = adims[1] = n;
289 :     SET_SLOT(ans, Matrix_uploSym, mkString("U"));
290 :     SET_SLOT(ans, Matrix_iSym,
291 :     Matrix_lgCsyrk(1, tra, n, k,
292 :     INTEGER(GET_SLOT(x, Matrix_iSym)),
293 :     INTEGER(GET_SLOT(x, Matrix_pSym)),
294 :     0, R_NilValue,
295 :     INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1))));
296 :     UNPROTECT(1);
297 :     return ans;
298 :     }

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