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 1064 - (view) (download) (as text)

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

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