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

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