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 956 - (view) (download) (as text)
Original Path: pkg/src/dgCMatrix.c

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

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