SCM

SCM Repository

[matrix] Annotation of /pkg/src/cscMatrix.c
ViewVC logotype

Annotation of /pkg/src/cscMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 310 - (view) (download) (as text)

1 : bates 10 #include "cscMatrix.h"
2 :    
3 :     SEXP csc_validate(SEXP x)
4 :     {
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 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sscMatrix"))), tmp;
45 :     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 :     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 :     return cscMatrix_set_Dim(ans, ncol);
101 :     }
102 :    
103 :     SEXP csc_matrix_crossprod(SEXP x, SEXP y)
104 :     {
105 :     SEXP pslot = GET_SLOT(x, Matrix_pSym), ans;
106 :     int j,
107 :     *xp = INTEGER(pslot),
108 :     *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
109 :     xncol = length(pslot) - 1,
110 :     xnrow = INTEGER(GET_SLOT(x, Matrix_DimSym))[0],
111 :     *ydims;
112 :     double *xx = REAL(GET_SLOT(x, Matrix_xSym));
113 :    
114 :     if (!(isMatrix(y) && isReal(y))) error("y must be a numeric matrix");
115 :     ydims = INTEGER(getAttrib(y, R_DimSymbol));
116 :     if (xnrow != ydims[0]) error("x and y must have the same number of rows");
117 :     ans = PROTECT(allocMatrix(REALSXP, xncol, ydims[1]));
118 :     for (j = 0; j < ydims[1]; j++) {
119 :     int i; double *ypt = REAL(y) + j * ydims[0];
120 :     for(i = 0; i < xncol; i++) {
121 :     int ii; double accum = 0.;
122 :     for (ii = xp[i]; ii < xp[i+1]; ii++) {
123 :     accum += xx[ii] * ypt[xi[ii]];
124 :     }
125 :     REAL(ans)[i + j * xncol] = accum;
126 :     }
127 :     }
128 :     UNPROTECT(1);
129 :     return ans;
130 :     }
131 :    
132 :     SEXP csc_to_triplet(SEXP x)
133 :     {
134 :     SEXP
135 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tripletMatrix"))),
136 :     dimslot = GET_SLOT(x, Matrix_DimSym),
137 :     islot = GET_SLOT(x, Matrix_iSym),
138 :     pslot = GET_SLOT(x, Matrix_pSym);
139 :     int *dims = INTEGER(dimslot), j, jj,
140 :     *xp = INTEGER(pslot), *yj;
141 :    
142 :     SET_SLOT(ans, Matrix_iSym, duplicate(islot));
143 :     SET_SLOT(ans, Matrix_DimSym, duplicate(dimslot));
144 :     SET_SLOT(ans, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));
145 :     SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, length(islot)));
146 :     yj = INTEGER(GET_SLOT(ans, Matrix_jSym));
147 :     jj = 0;
148 :     for (j = 0; j < dims[1]; j++) {
149 :     while (jj < xp[j + 1]) {
150 :     yj[jj] = j;
151 :     jj++;
152 :     }
153 :     }
154 :     UNPROTECT(1);
155 :     return ans;
156 :     }
157 :    
158 :     SEXP csc_to_matrix(SEXP x)
159 :     {
160 :     SEXP ans, pslot = GET_SLOT(x, Matrix_pSym);
161 :     int j, ncol = length(pslot) - 1,
162 :     nrow = INTEGER(GET_SLOT(x, Matrix_DimSym))[0],
163 :     *xp = INTEGER(pslot),
164 :     *xi = INTEGER(GET_SLOT(x, Matrix_iSym));
165 :     double *xx = REAL(GET_SLOT(x, Matrix_xSym)), *ax;
166 :    
167 :     ax = REAL(ans = PROTECT(allocMatrix(REALSXP, nrow, ncol)));
168 :     for (j = 0; j < (nrow * ncol); j++) ax[j] = 0.;
169 :     for (j = 0; j < ncol; j++) {
170 :     int ind;
171 :     for (ind = xp[j]; ind < xp[j+1]; ind++) {
172 :     ax[j * nrow + xi[ind]] = xx[ind];
173 :     }
174 :     }
175 :     UNPROTECT(1);
176 :     return ans;
177 :     }
178 :    
179 :     SEXP csc_to_geMatrix(SEXP x)
180 :     {
181 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("geMatrix"))),
182 :     Dimslot = GET_SLOT(x, Matrix_DimSym);
183 :     int *dims = INTEGER(Dimslot),
184 :     *xp = INTEGER(GET_SLOT(x, Matrix_pSym)),
185 :     *xi = INTEGER(GET_SLOT(x, Matrix_iSym));
186 :     double *xx = REAL(GET_SLOT(x, Matrix_xSym)), *ax;
187 :     int j, nrow = dims[0], ncol = dims[1];
188 :    
189 :     SET_SLOT(ans, Matrix_DimSym, duplicate(Dimslot));
190 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nrow*ncol));
191 :     ax = REAL(GET_SLOT(ans, Matrix_xSym));
192 :     for (j = 0; j < (nrow * ncol); j++) ax[j] = 0.;
193 :     for (j = 0; j < ncol; j++) {
194 :     int ind;
195 :     for (ind = xp[j]; ind < xp[j+1]; ind++) {
196 :     ax[j * nrow + xi[ind]] = xx[ind];
197 :     }
198 :     }
199 :     UNPROTECT(1);
200 :     return ans;
201 :     }
202 :    
203 : bates 310 /* FIXME: This function does not appear to be used. */
204 :     static /* added to check if it was used */
205 : bates 10 SEXP csc_to_imagemat(SEXP x)
206 :     {
207 :     SEXP ans, pslot = GET_SLOT(x, Matrix_pSym);
208 :     int j, ncol = length(pslot) - 1,
209 :     nrow = INTEGER(GET_SLOT(x, Matrix_DimSym))[0],
210 :     *xp = INTEGER(pslot),
211 :     *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
212 :     *ax;
213 :     /* ans is the transpose of the indicator
214 :     of non-zero */
215 :     ax = INTEGER(ans = PROTECT(allocMatrix(INTSXP, ncol, nrow)));
216 :     for (j = 0; j < (ncol * nrow); j++) ax[j] = 0;
217 :     for (j = 0; j < ncol; j++) {
218 :     int ind;
219 :    
220 :     for (ind = xp[j]; ind < xp[j+1]; ind++) {
221 :     /* reverse rows of transpose */
222 :     ax[j + ncol*(nrow - 1 - xi[ind])] = 1;
223 :     }
224 :     }
225 :     UNPROTECT(1);
226 :     return ans;
227 :     }
228 :    
229 :     SEXP matrix_to_csc(SEXP A)
230 :     {
231 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("cscMatrix")));
232 :     int *adims = INTEGER(getAttrib(A, R_DimSymbol)), j,
233 :     maxnz, nrow, ncol, nnz, *vp, *vi;
234 :    
235 :     double *vx;
236 :    
237 :     if (!(isMatrix(A) && isReal(A)))
238 :     error("A must be a numeric matrix");
239 :     nrow = adims[0]; ncol = adims[1];
240 :     SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, ncol + 1));
241 :     vp = INTEGER(GET_SLOT(val, Matrix_pSym));
242 :     maxnz = nrow * ncol;
243 :     vi = Calloc(maxnz, int); vx = Calloc(maxnz, double);
244 :     nnz = 0;
245 :     for (j = 0; j < ncol; j++) {
246 :     int i;
247 :     vp[j] = nnz;
248 :     for (i = 0; i < nrow; i++) {
249 :     double val = REAL(A)[i + j * nrow];
250 :     if (val != 0.) {
251 :     vi[nnz] = i;
252 :     vx[nnz] = val;
253 :     nnz++;
254 :     }
255 :     }
256 :     }
257 :     vp[ncol] = nnz;
258 :     SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));
259 :     Memcpy(INTEGER(GET_SLOT(val, Matrix_iSym)), vi, nnz);
260 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));
261 :     Memcpy(REAL(GET_SLOT(val, Matrix_xSym)), vx, nnz);
262 :     Free(vi); Free(vx);
263 :     UNPROTECT(1);
264 :     return cscMatrix_set_Dim(val, nrow);
265 :     }
266 :    
267 :    
268 :     SEXP triplet_to_csc(SEXP triplet)
269 :     {
270 :     SEXP Tisl = GET_SLOT(triplet, Matrix_iSym);
271 :     int *Ti = INTEGER(Tisl),
272 :     *Tj = INTEGER(GET_SLOT(triplet, Matrix_jSym)),
273 :     i, nrow, ncol,
274 :     nz = length(Tisl);
275 :    
276 :     nrow = ncol = -1;
277 :     for(i = 0; i < nz; i++) {
278 :     if (Ti[i] > nrow) nrow = Ti[i];
279 :     if (Tj[i] > ncol) ncol = Tj[i];
280 :     }
281 :     return triple_as_SEXP(nrow + 1, ncol + 1, nz, Ti, Tj,
282 :     REAL(GET_SLOT(triplet, Matrix_xSym)),
283 :     "cscMatrix");
284 :     }
285 :    
286 :     SEXP csc_getDiag(SEXP x)
287 :     {
288 :     SEXP pslot = GET_SLOT(x, Matrix_pSym), ans;
289 :     int *xp = INTEGER(pslot),
290 :     *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
291 :     j,
292 :     ncol = length(pslot) - 1,
293 :     nrow = INTEGER(GET_SLOT(x, Matrix_DimSym))[0],
294 :     ndiag;
295 :     double *xx = REAL(GET_SLOT(x, Matrix_xSym)), *diag;
296 :    
297 :     ndiag = (nrow < ncol) ? nrow : ncol;
298 :     ans = PROTECT(allocVector(REALSXP, ndiag));
299 :     diag = REAL(ans);
300 :     for (j = 0; j < ndiag; j++) {
301 :     int ind;
302 :     diag[j] = 0.;
303 :     for (ind = xp[j]; ind < xp[j+1]; ind++) {
304 :     if (xi[ind] == j) diag[j] = xx[ind];
305 :     }
306 :     }
307 :     UNPROTECT(1);
308 :     return ans;
309 :     }
310 :    
311 :     SEXP csc_transpose(SEXP x)
312 :     {
313 :     SEXP
314 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS("cscMatrix"))),
315 :     islot = GET_SLOT(x, Matrix_iSym);
316 :     int nnz = length(islot),
317 :     *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),
318 :     *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
319 :    
320 :     adims[0] = xdims[1]; adims[1] = xdims[0];
321 :     SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, xdims[0] + 1));
322 :     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
323 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
324 :     csc_components_transpose(xdims[0], xdims[1], nnz,
325 :     INTEGER(GET_SLOT(x, Matrix_pSym)),
326 :     INTEGER(islot),
327 :     REAL(GET_SLOT(x, Matrix_xSym)),
328 :     INTEGER(GET_SLOT(ans, Matrix_pSym)),
329 :     INTEGER(GET_SLOT(ans, Matrix_iSym)),
330 :     REAL(GET_SLOT(ans, Matrix_xSym)));
331 :     UNPROTECT(1);
332 :     return ans;
333 :     }
334 : bates 310
335 :     SEXP csc_matrix_mm(SEXP a, SEXP b)
336 :     {
337 :     int *adim = INTEGER(GET_SLOT(a, Matrix_DimSym)),
338 :     *ai = INTEGER(GET_SLOT(a, Matrix_iSym)),
339 :     *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),
340 :     *bdim = INTEGER(getAttrib(b, R_DimSymbol));
341 :     int j, k, m = adim[0], n = bdim[1], r = adim[1];
342 :     double *ax = REAL(GET_SLOT(a, Matrix_xSym));
343 :     SEXP val;
344 :    
345 :     if (bdim[0] != r)
346 :     error("Matrices of sizes (%d,%d) and (%d,%d) cannot be multiplied",
347 :     m, r, bdim[0], n);
348 :     val = PROTECT(allocMatrix(REALSXP, m, n));
349 :     for (j = 0; j < n; j++) { /* across columns of b */
350 :     double *ccol = REAL(val) + j * m,
351 :     *bcol = REAL(b) + j * r;
352 :    
353 :     for (k = 0; k < m; k++) ccol[k] = 0.; /* zero the accumulators */
354 :     for (k = 0; k < r; k++) { /* across columns of a */
355 :     int kk, k2 = ap[k + 1];
356 :     for (kk = ap[k]; kk < k2; kk++) {
357 :     ccol[ai[kk]] += ax[kk] * bcol[k];
358 :     }
359 :     }
360 :     }
361 :     UNPROTECT(1);
362 :     return val;
363 :     }
364 :    
365 :    
366 :    

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