SCM

SCM Repository

[matrix] Annotation of /branches/Matrix-mer2/src/dgCMatrix.c
ViewVC logotype

Annotation of /branches/Matrix-mer2/src/dgCMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 677 - (view) (download) (as text)
Original Path: pkg/src/dgCMatrix.c

1 : bates 478 #include "dgCMatrix.h"
2 : bates 10
3 : bates 488 SEXP dgCMatrix_validate(SEXP x)
4 : bates 10 {
5 :     SEXP pslot = GET_SLOT(x, Matrix_pSym),
6 :     islot = GET_SLOT(x, Matrix_iSym),
7 :     xslot = GET_SLOT(x, Matrix_xSym);
8 :     int j,
9 :     ncol = length(pslot) - 1,
10 :     *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
11 :     nrow,
12 :     *xp = INTEGER(pslot),
13 :     *xi = INTEGER(islot);
14 : maechler 534
15 : bates 10 nrow = dims[0];
16 :     if (length(islot) != length(xslot))
17 : bates 582 return mkString(_("lengths of slots i and x must match"));
18 : bates 10 if (length(pslot) <= 0)
19 : bates 582 return mkString(_("slot p must have length > 0"));
20 : bates 10 if (xp[0] != 0)
21 : bates 582 return mkString(_("first element of slot p must be zero"));
22 : bates 10 if (length(islot) != xp[ncol])
23 : bates 582 return mkString(_("last element of slot p must match length of slots i and x"));
24 : bates 10 for (j = 0; j < ncol; j++) {
25 :     if (xp[j] > xp[j+1])
26 : bates 582 return mkString(_("slot p must be non-decreasing"));
27 : bates 10 }
28 :     for (j = 0; j < length(islot); j++) {
29 :     if (xi[j] < 0 || xi[j] >= nrow)
30 : bates 582 return mkString(_("all row indices must be between 0 and nrow-1"));
31 : bates 10 }
32 :     if (csc_unsorted_columns(ncol, xp, xi)) {
33 :     csc_sort_columns(ncol, xp, xi, REAL(xslot));
34 :     }
35 :     return ScalarLogical(1);
36 :     }
37 : maechler 534
38 : bates 10 SEXP csc_crossprod(SEXP x)
39 :     {
40 :     SEXP pslot = GET_SLOT(x, Matrix_pSym),
41 : bates 478 ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsCMatrix"))), tmp;
42 : bates 10 int *xp = INTEGER(pslot),
43 :     *xi = INTEGER(GET_SLOT(x, Matrix_iSym));
44 :     double *xx = REAL(GET_SLOT(x, Matrix_xSym));
45 :    
46 :     int j, *iVal, ncol = length(pslot) - 1, maxnz, nnz = 0, *pVal;
47 :     double *xVal;
48 : maechler 534
49 : bates 476 SET_SLOT(ans, Matrix_factorSym, allocVector(VECSXP, 0));
50 : bates 492 SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
51 :     SET_SLOT(ans, Matrix_uploSym, mkString("L"));
52 : bates 10 maxnz = (ncol * (ncol + 1))/2;
53 :     iVal = Calloc(maxnz, int); xVal = Calloc(maxnz, double);
54 :     SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, ncol + 1));
55 :     tmp = GET_SLOT(ans, Matrix_pSym);
56 :     pVal = INTEGER(tmp);
57 :     for (j = 0; j < ncol; j++) {
58 :     pVal[j] = nnz;
59 :     if (xp[j] < xp[j+1]) { /* column j contains some non-zeros */
60 :     int ind, jj;
61 :     double accum = 0.;
62 :     /* diagonal elements */
63 :     for (ind = xp[j]; ind < xp[j+1]; ind++)
64 :     accum += xx[ind] * xx[ind];
65 :     iVal[nnz] = j;
66 :     xVal[nnz] = accum;
67 :     nnz++;
68 :     /* off-diagonals (lower triangle only) */
69 :     for (jj = j+1; jj < ncol; jj++) {
70 :     int ind2;
71 :    
72 :     ind = xp[j];
73 :     ind2 = xp[jj];
74 :     accum = 0.;
75 :     while (ind < xp[j+1] && ind2 < xp[jj+1]) {
76 :     if (xi[ind] < xi[ind2]) ind++;
77 :     else {
78 :     if (xi[ind] > xi[ind2]) ind2++;
79 :     else {
80 :     accum += xx[ind] * xx[ind2];
81 :     ind++; ind2++;
82 :     }
83 :     }
84 :     }
85 :     if (accum != 0.) {
86 :     iVal[nnz] = jj;
87 :     xVal[nnz] = accum;
88 :     nnz++;
89 :     }
90 :     }
91 :     }
92 :     }
93 :     pVal[ncol] = nnz;
94 :    
95 :     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
96 :     Memcpy(INTEGER(GET_SLOT(ans, Matrix_iSym)), iVal, nnz);
97 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
98 :     Memcpy(REAL(GET_SLOT(ans, Matrix_xSym)), xVal, nnz);
99 :     Free(iVal); Free(xVal); UNPROTECT(1);
100 : bates 478 return dgCMatrix_set_Dim(ans, ncol);
101 : bates 10 }
102 :    
103 : bates 332 SEXP csc_tcrossprod(SEXP x)
104 :     {
105 :     SEXP pslot = GET_SLOT(x, Matrix_pSym),
106 : bates 478 ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsCMatrix")));
107 : bates 332 int *xp = INTEGER(pslot),
108 :     *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
109 :     *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
110 :     double *xx = REAL(GET_SLOT(x, Matrix_xSym));
111 :    
112 :     int j, ntrip, *iVal, nrow = dims[0], ncol = dims[1], *jVal, nnz, pos;
113 :     int *itmp, *ansp;
114 :     double *xVal, *xtmp;
115 :    
116 : bates 476 SET_SLOT(ans, Matrix_factorSym, allocVector(VECSXP, 0));
117 : bates 492 SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
118 : bates 332 ntrip = nrow; /* number of triplets */
119 :     for (j = 0; j < ncol; j++) {
120 :     int nzj = xp[j+1] - xp[j];
121 :     ntrip += (nzj * (nzj - 1))/2;
122 :     }
123 :     iVal = Calloc(ntrip, int); jVal = Calloc(ntrip, int);
124 :     xVal = Calloc(ntrip, double);
125 :     for (j = 0; j < nrow; j++) {
126 :     iVal[j] = jVal[j] = j;
127 :     xVal[j] = 0.;
128 :     }
129 :     pos = nrow;
130 :     for (j = 0; j < ncol; j++) {
131 :     int k, kk, k2 = xp[j+1];
132 :     for (k = xp[j]; k < k2; k++) {
133 :     int r1 = xi[k];
134 :     double x1 = xx[k];
135 :     xVal[r1] += x1 * x1;
136 :     for (kk = k + 1; kk < k2; kk++) {
137 :     int r2 = xi[kk];
138 :     double x2 = xx[kk];
139 :     jVal[pos] = r1;
140 :     iVal[pos] = r2;
141 :     xVal[pos] = x1 * x2;
142 :     pos++;
143 :     }
144 :     }
145 :     }
146 :     SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, nrow + 1));
147 :     ansp = INTEGER(GET_SLOT(ans, Matrix_pSym));
148 :     itmp = Calloc(ntrip, int); xtmp = Calloc(ntrip, double);
149 : bates 488 triplet_to_col(nrow, nrow, ntrip, iVal, jVal, xVal,
150 : bates 332 ansp, itmp, xtmp);
151 :     nnz = ansp[nrow];
152 : maechler 534 SET_SLOT(ans, Matrix_uploSym, mkString("L"));
153 : bates 332 SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
154 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
155 :     Memcpy(INTEGER(GET_SLOT(ans, Matrix_iSym)), itmp, nnz);
156 :     Memcpy(REAL(GET_SLOT(ans, Matrix_xSym)), xtmp, nnz);
157 :     dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
158 :     dims[0] = dims[1] = nrow;
159 :     Free(itmp); Free(xtmp); Free(iVal); Free(jVal); Free(xVal);
160 :     UNPROTECT(1);
161 :     return ans;
162 :     }
163 : maechler 534
164 : bates 10 SEXP csc_matrix_crossprod(SEXP x, SEXP y)
165 :     {
166 :     SEXP pslot = GET_SLOT(x, Matrix_pSym), ans;
167 :     int j,
168 :     *xp = INTEGER(pslot),
169 :     *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
170 :     xncol = length(pslot) - 1,
171 :     xnrow = INTEGER(GET_SLOT(x, Matrix_DimSym))[0],
172 :     *ydims;
173 :     double *xx = REAL(GET_SLOT(x, Matrix_xSym));
174 :    
175 : bates 582 if (!(isMatrix(y) && isReal(y))) error(_("y must be a numeric matrix"));
176 : bates 10 ydims = INTEGER(getAttrib(y, R_DimSymbol));
177 : bates 582 if (xnrow != ydims[0]) error(_("x and y must have the same number of rows"));
178 : bates 10 ans = PROTECT(allocMatrix(REALSXP, xncol, ydims[1]));
179 :     for (j = 0; j < ydims[1]; j++) {
180 :     int i; double *ypt = REAL(y) + j * ydims[0];
181 :     for(i = 0; i < xncol; i++) {
182 :     int ii; double accum = 0.;
183 :     for (ii = xp[i]; ii < xp[i+1]; ii++) {
184 :     accum += xx[ii] * ypt[xi[ii]];
185 :     }
186 :     REAL(ans)[i + j * xncol] = accum;
187 :     }
188 :     }
189 :     UNPROTECT(1);
190 :     return ans;
191 :     }
192 :    
193 : bates 677 SEXP compressed_to_dgTMatrix(SEXP x, SEXP colP)
194 : bates 10 {
195 : bates 677 int col = asLogical(colP);
196 :     SEXP indSym = col ? Matrix_iSym : Matrix_jSym;
197 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix"))),
198 :     indP = GET_SLOT(x, indSym),
199 :     pP = GET_SLOT(x, Matrix_pSym);
200 :     int npt = length(pP);
201 : maechler 534
202 : bates 677 SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
203 : bates 10 SET_SLOT(ans, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));
204 : bates 677 SET_SLOT(ans, indSym, duplicate(indP));
205 :     expand_column_pointers(length(pP) - 1, INTEGER(pP),
206 :     INTEGER(ALLOC_SLOT(ans,
207 :     col ? Matrix_jSym : Matrix_iSym,
208 :     INTSXP, length(indP))));
209 : bates 10 UNPROTECT(1);
210 :     return ans;
211 :     }
212 :    
213 :     SEXP csc_to_matrix(SEXP x)
214 :     {
215 :     SEXP ans, pslot = GET_SLOT(x, Matrix_pSym);
216 :     int j, ncol = length(pslot) - 1,
217 :     nrow = INTEGER(GET_SLOT(x, Matrix_DimSym))[0],
218 :     *xp = INTEGER(pslot),
219 :     *xi = INTEGER(GET_SLOT(x, Matrix_iSym));
220 :     double *xx = REAL(GET_SLOT(x, Matrix_xSym)), *ax;
221 :    
222 :     ax = REAL(ans = PROTECT(allocMatrix(REALSXP, nrow, ncol)));
223 :     for (j = 0; j < (nrow * ncol); j++) ax[j] = 0.;
224 :     for (j = 0; j < ncol; j++) {
225 :     int ind;
226 :     for (ind = xp[j]; ind < xp[j+1]; ind++) {
227 :     ax[j * nrow + xi[ind]] = xx[ind];
228 :     }
229 :     }
230 :     UNPROTECT(1);
231 :     return ans;
232 :     }
233 :    
234 : bates 478 SEXP csc_to_dgeMatrix(SEXP x)
235 : bates 10 {
236 : bates 478 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
237 : bates 10 Dimslot = GET_SLOT(x, Matrix_DimSym);
238 :     int *dims = INTEGER(Dimslot),
239 :     *xp = INTEGER(GET_SLOT(x, Matrix_pSym)),
240 :     *xi = INTEGER(GET_SLOT(x, Matrix_iSym));
241 :     double *xx = REAL(GET_SLOT(x, Matrix_xSym)), *ax;
242 :     int j, nrow = dims[0], ncol = dims[1];
243 : maechler 534
244 : bates 10 SET_SLOT(ans, Matrix_DimSym, duplicate(Dimslot));
245 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nrow*ncol));
246 : bates 332 SET_SLOT(ans, Matrix_rcondSym, allocVector(REALSXP, 0));
247 : bates 476 SET_SLOT(ans, Matrix_factorSym, allocVector(VECSXP, 0));
248 : bates 10 ax = REAL(GET_SLOT(ans, Matrix_xSym));
249 :     for (j = 0; j < (nrow * ncol); j++) ax[j] = 0.;
250 :     for (j = 0; j < ncol; j++) {
251 :     int ind;
252 :     for (ind = xp[j]; ind < xp[j+1]; ind++) {
253 :     ax[j * nrow + xi[ind]] = xx[ind];
254 :     }
255 :     }
256 :     UNPROTECT(1);
257 :     return ans;
258 :     }
259 :    
260 :     SEXP matrix_to_csc(SEXP A)
261 :     {
262 : bates 478 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix")));
263 : bates 10 int *adims = INTEGER(getAttrib(A, R_DimSymbol)), j,
264 :     maxnz, nrow, ncol, nnz, *vp, *vi;
265 : maechler 534
266 : bates 10 double *vx;
267 :    
268 :     if (!(isMatrix(A) && isReal(A)))
269 : bates 582 error(_("A must be a numeric matrix"));
270 : bates 10 nrow = adims[0]; ncol = adims[1];
271 : bates 476 SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
272 : bates 492 SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
273 : bates 10 SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, ncol + 1));
274 :     vp = INTEGER(GET_SLOT(val, Matrix_pSym));
275 :     maxnz = nrow * ncol;
276 :     vi = Calloc(maxnz, int); vx = Calloc(maxnz, double);
277 :     nnz = 0;
278 :     for (j = 0; j < ncol; j++) {
279 :     int i;
280 :     vp[j] = nnz;
281 :     for (i = 0; i < nrow; i++) {
282 :     double val = REAL(A)[i + j * nrow];
283 :     if (val != 0.) {
284 :     vi[nnz] = i;
285 :     vx[nnz] = val;
286 :     nnz++;
287 :     }
288 :     }
289 :     }
290 :     vp[ncol] = nnz;
291 :     SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));
292 :     Memcpy(INTEGER(GET_SLOT(val, Matrix_iSym)), vi, nnz);
293 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));
294 :     Memcpy(REAL(GET_SLOT(val, Matrix_xSym)), vx, nnz);
295 :     Free(vi); Free(vx);
296 :     UNPROTECT(1);
297 : bates 478 return dgCMatrix_set_Dim(val, nrow);
298 : bates 10 }
299 :    
300 : maechler 534
301 : bates 478 SEXP dgTMatrix_to_csc(SEXP dgTMatrix)
302 : bates 10 {
303 : bates 478 SEXP Tisl = GET_SLOT(dgTMatrix, Matrix_iSym);
304 : bates 10 int *Ti = INTEGER(Tisl),
305 : bates 478 *Tj = INTEGER(GET_SLOT(dgTMatrix, Matrix_jSym)),
306 : bates 10 i, nrow, ncol,
307 :     nz = length(Tisl);
308 :    
309 :     nrow = ncol = -1;
310 :     for(i = 0; i < nz; i++) {
311 :     if (Ti[i] > nrow) nrow = Ti[i];
312 :     if (Tj[i] > ncol) ncol = Tj[i];
313 :     }
314 :     return triple_as_SEXP(nrow + 1, ncol + 1, nz, Ti, Tj,
315 : bates 478 REAL(GET_SLOT(dgTMatrix, Matrix_xSym)),
316 :     "dgCMatrix");
317 : bates 10 }
318 :    
319 :     SEXP csc_getDiag(SEXP x)
320 :     {
321 :     SEXP pslot = GET_SLOT(x, Matrix_pSym), ans;
322 :     int *xp = INTEGER(pslot),
323 :     *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
324 :     j,
325 :     ncol = length(pslot) - 1,
326 :     nrow = INTEGER(GET_SLOT(x, Matrix_DimSym))[0],
327 :     ndiag;
328 :     double *xx = REAL(GET_SLOT(x, Matrix_xSym)), *diag;
329 :    
330 :     ndiag = (nrow < ncol) ? nrow : ncol;
331 :     ans = PROTECT(allocVector(REALSXP, ndiag));
332 :     diag = REAL(ans);
333 :     for (j = 0; j < ndiag; j++) {
334 :     int ind;
335 :     diag[j] = 0.;
336 :     for (ind = xp[j]; ind < xp[j+1]; ind++) {
337 :     if (xi[ind] == j) diag[j] = xx[ind];
338 :     }
339 :     }
340 :     UNPROTECT(1);
341 :     return ans;
342 :     }
343 :    
344 :     SEXP csc_transpose(SEXP x)
345 :     {
346 : bates 492 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix"))),
347 : bates 10 islot = GET_SLOT(x, Matrix_iSym);
348 : bates 492 int *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
349 :     nnz = length(islot);
350 : bates 10
351 : bates 588 adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
352 : bates 10 adims[0] = xdims[1]; adims[1] = xdims[0];
353 : bates 588 csc_compTr(xdims[0], xdims[1], nnz,
354 :     INTEGER(GET_SLOT(x, Matrix_pSym)), INTEGER(islot),
355 :     REAL(GET_SLOT(x, Matrix_xSym)),
356 :     INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, xdims[0] + 1)),
357 :     INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)),
358 :     REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)));
359 : bates 10 UNPROTECT(1);
360 :     return ans;
361 :     }
362 : bates 310
363 :     SEXP csc_matrix_mm(SEXP a, SEXP b)
364 :     {
365 :     int *adim = INTEGER(GET_SLOT(a, Matrix_DimSym)),
366 :     *ai = INTEGER(GET_SLOT(a, Matrix_iSym)),
367 :     *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),
368 :     *bdim = INTEGER(getAttrib(b, R_DimSymbol));
369 :     int j, k, m = adim[0], n = bdim[1], r = adim[1];
370 :     double *ax = REAL(GET_SLOT(a, Matrix_xSym));
371 :     SEXP val;
372 :    
373 :     if (bdim[0] != r)
374 : bates 582 error(_("Matrices of sizes (%d,%d) and (%d,%d) cannot be multiplied"),
375 : bates 310 m, r, bdim[0], n);
376 :     val = PROTECT(allocMatrix(REALSXP, m, n));
377 :     for (j = 0; j < n; j++) { /* across columns of b */
378 :     double *ccol = REAL(val) + j * m,
379 :     *bcol = REAL(b) + j * r;
380 :    
381 :     for (k = 0; k < m; k++) ccol[k] = 0.; /* zero the accumulators */
382 :     for (k = 0; k < r; k++) { /* across columns of a */
383 :     int kk, k2 = ap[k + 1];
384 :     for (kk = ap[k]; kk < k2; kk++) {
385 :     ccol[ai[kk]] += ax[kk] * bcol[k];
386 :     }
387 :     }
388 :     }
389 :     UNPROTECT(1);
390 :     return val;
391 :     }
392 :    
393 : bates 332 SEXP csc_col_permute(SEXP x, SEXP perm)
394 :     {
395 : bates 478 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix"))), tmp;
396 : bates 332 int *iperm, *prm, *vi, *vp, *xi, *xp, j, k, ncol, pos;
397 :     double *vx, *xx;
398 : bates 310
399 : bates 476 SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
400 : bates 332 tmp = GET_SLOT(x, Matrix_DimSym);
401 :     SET_SLOT(val, Matrix_DimSym, duplicate(tmp));
402 :     ncol = INTEGER(tmp)[1];
403 :     if (!(isInteger(perm) && length(perm) == ncol))
404 : bates 582 error(_("perm must be an integer vector of length %d"),
405 : bates 332 ncol);
406 :     prm = INTEGER(perm);
407 : bates 371 if (!R_ldl_valid_perm(ncol, prm))
408 : bates 582 error(_("perm is not a valid 0-based permutation"));
409 : bates 332 iperm = Calloc(ncol, int);
410 :     for (j = 0; j < ncol; j++) iperm[prm[j]] = j;
411 :     tmp = GET_SLOT(x, Matrix_pSym);
412 :     xp = INTEGER(tmp);
413 :     SET_SLOT(val, Matrix_pSym, duplicate(tmp));
414 :     vp = INTEGER(GET_SLOT(val, Matrix_pSym));
415 :     tmp = GET_SLOT(x, Matrix_iSym);
416 :     xi = INTEGER(tmp);
417 :     SET_SLOT(val, Matrix_iSym, duplicate(tmp));
418 :     vi = INTEGER(GET_SLOT(val, Matrix_iSym));
419 :     tmp = GET_SLOT(x, Matrix_xSym);
420 :     xx = REAL(tmp);
421 :     SET_SLOT(val, Matrix_xSym, duplicate(tmp));
422 :     vx = REAL(GET_SLOT(val, Matrix_xSym));
423 :    
424 :     pos = vp[0] = 0;
425 :     for (j = 0; j < ncol; j++) {
426 :     int jj = iperm[j];
427 :     int j1 = xp[jj], j2 = xp[jj+1];
428 :     vp[j + 1] = vp[j] + (j2 - j1);
429 :     for (k = j1; k < j2; k++) {
430 :     vi[pos] = xi[k];
431 :     vx[pos] = xx[k];
432 :     pos++;
433 :     }
434 :     }
435 :     Free(iperm);
436 :     UNPROTECT(1);
437 :     return val;
438 :     }
439 :    
440 : maechler 534
441 :    

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