SCM

SCM Repository

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

Annotation of /pkg/Matrix/src/dtCMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : bates 741 /* Sparse triangular numeric matrices */
2 : bates 478 #include "dtCMatrix.h"
3 : bates 1251 #include "cs_utils.h"
4 : bates 10
5 :     SEXP tsc_validate(SEXP x)
6 :     {
7 : maechler 890 return triangularMatrix_validate(x);
8 : maechler 895 /* see ./dsCMatrix.c or ./dtpMatrix.c on how to do more testing here */
9 : bates 10 }
10 :    
11 : bates 1209 #if 0
12 : bates 10 SEXP tsc_transpose(SEXP x)
13 :     {
14 : bates 1209 cholmod_sparse *cx = as_cholmod_sparse(x);
15 :    
16 : maechler 945 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix"))),
17 : bates 10 islot = GET_SLOT(x, Matrix_iSym);
18 :     int nnz = length(islot),
19 : bates 338 *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
20 : maechler 951 int up = uplo_P(x)[0] == 'U';
21 : bates 10
22 : bates 588 adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
23 : bates 10 adims[0] = xdims[1]; adims[1] = xdims[0];
24 : maechler 945
25 : maechler 951 if(*diag_P(x) == 'U')
26 : maechler 945 SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(x, Matrix_diagSym)));
27 : bates 741 SET_SLOT(ans, Matrix_uploSym, mkString(up ? "L" : "U"));
28 : maechler 945
29 : bates 588 csc_compTr(xdims[0], xdims[1], nnz,
30 :     INTEGER(GET_SLOT(x, Matrix_pSym)), INTEGER(islot),
31 :     REAL(GET_SLOT(x, Matrix_xSym)),
32 :     INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, xdims[0] + 1)),
33 :     INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)),
34 :     REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)));
35 : bates 10 UNPROTECT(1);
36 :     return ans;
37 :     }
38 : bates 1209 #endif
39 : bates 70
40 : bates 478 SEXP tsc_to_dgTMatrix(SEXP x)
41 : bates 70 {
42 :     SEXP ans;
43 : maechler 951 if (*diag_P(x) != 'U')
44 : bates 677 ans = compressed_to_dgTMatrix(x, ScalarLogical(1));
45 : bates 70 else { /* unit triangular matrix */
46 : maechler 534 SEXP islot = GET_SLOT(x, Matrix_iSym),
47 : bates 70 pslot = GET_SLOT(x, Matrix_pSym);
48 :     int *ai, *aj, j,
49 :     n = length(pslot) - 1,
50 :     nod = length(islot),
51 :     nout = n + nod,
52 :     *p = INTEGER(pslot);
53 :     double *ax;
54 : maechler 534
55 : bates 478 ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix")));
56 : bates 70 SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
57 :     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nout));
58 :     ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
59 :     Memcpy(ai, INTEGER(islot), nod);
60 :     SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, nout));
61 :     aj = INTEGER(GET_SLOT(ans, Matrix_jSym));
62 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nout));
63 :     ax = REAL(GET_SLOT(ans, Matrix_xSym));
64 :     Memcpy(ax, REAL(GET_SLOT(x, Matrix_xSym)), nod);
65 :     for (j = 0; j < n; j++) {
66 :     int jj, npj = nod + j, p2 = p[j+1];
67 :     aj[npj] = ai[npj] = j;
68 :     ax[npj] = 1.;
69 :     for (jj = p[j]; jj < p2; jj++) aj[jj] = j;
70 :     }
71 :     UNPROTECT(1);
72 :     }
73 :     return ans;
74 :     }
75 : bates 338
76 : maechler 534 /**
77 : bates 358 * Derive the column pointer vector for the inverse of L from the parent array
78 : maechler 534 *
79 : bates 358 * @param n length of parent array
80 :     * @param countDiag 0 for a unit triangular matrix with implicit diagonal, otherwise 1
81 :     * @param pr parent vector describing the elimination tree
82 :     * @param ap array of length n+1 to be filled with the column pointers
83 : maechler 534 *
84 : bates 358 * @return the number of non-zero entries (ap[n])
85 :     */
86 :     int parent_inv_ap(int n, int countDiag, const int pr[], int ap[])
87 :     {
88 :     int *sz = Calloc(n, int), j;
89 :    
90 :     for (j = n - 1; j >= 0; j--) {
91 :     int parent = pr[j];
92 :     sz[j] = (parent < 0) ? countDiag : (1 + sz[parent]);
93 :     }
94 :     ap[0] = 0;
95 :     for (j = 0; j < n; j++)
96 :     ap[j+1] = ap[j] + sz[j];
97 :     Free(sz);
98 :     return ap[n];
99 :     }
100 :    
101 : maechler 534 /**
102 : bates 358 * Derive the row index array for the inverse of L from the parent array
103 : maechler 534 *
104 : bates 358 * @param n length of parent array
105 :     * @param countDiag 0 for a unit triangular matrix with implicit diagonal, otherwise 1
106 :     * @param pr parent vector describing the elimination tree
107 :     * @param ai row index vector of length ap[n]
108 :     */
109 :     void parent_inv_ai(int n, int countDiag, const int pr[], int ai[])
110 :     {
111 :     int j, k, pos = 0;
112 :     for (j = 0; j < n; j++) {
113 :     if (countDiag) ai[pos++] = j;
114 :     for (k = pr[j]; k >= 0; k = pr[k]) ai[pos++] = k;
115 :     }
116 :     }
117 : maechler 534
118 : bates 338 SEXP Parent_inverse(SEXP par, SEXP unitdiag)
119 :     {
120 : bates 478 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
121 : bates 338 int *ap, *ai, *dims, *pr = INTEGER(par),
122 :     countDiag = 1 - asLogical(unitdiag),
123 : bates 367 j, n = length(par), nnz;
124 : bates 338 double *ax;
125 : maechler 534
126 : bates 582 if (!isInteger(par)) error(_("par argument must be an integer vector"));
127 : bates 338 SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, n + 1));
128 :     ap = INTEGER(GET_SLOT(ans, Matrix_pSym));
129 : bates 358 nnz = parent_inv_ap(n, countDiag, pr, ap);
130 : bates 338 SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
131 :     ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
132 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
133 :     ax = REAL(GET_SLOT(ans, Matrix_xSym));
134 :     for (j = 0; j < nnz; j++) ax[j] = 1.;
135 :     SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
136 :     dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
137 :     dims[0] = dims[1] = n;
138 : maechler 534 SET_SLOT(ans, Matrix_uploSym, mkString("L"));
139 :     SET_SLOT(ans, Matrix_diagSym, (countDiag ? mkString("N") : mkString("U")));
140 : bates 358 parent_inv_ai(n, countDiag, pr, ai);
141 : bates 338 UNPROTECT(1);
142 :     return ans;
143 :     }
144 : bates 1248
145 : bates 1251 #if 0
146 : bates 1248 SEXP dtCMatrix_solve(SEXP a)
147 :     {
148 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
149 :     int lo = uplo_P(a)[0] == 'L', unit = diag_P(a)[0] == 'U',
150 :     n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0], nnz,
151 :     *ai = INTEGER(GET_SLOT(a,Matrix_iSym)),
152 :     *ap = INTEGER(GET_SLOT(a, Matrix_pSym)), *bi,
153 :     *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
154 :     int bnz = 10 * ap[n]; /* initial estimate of nnz in b */
155 :     int *ri = Calloc(bnz, int), *ind = Calloc(n, int), j;
156 :     double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *bx;
157 :    
158 :     SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
159 :     SET_SLOT(ans, Matrix_DimNamesSym,
160 :     duplicate(GET_SLOT(a, Matrix_DimNamesSym)));
161 :     SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
162 :     SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
163 :    
164 :     if (!(lo && unit))
165 :     error("code for non-unit or upper triangular not yet written");
166 :     /* Initially bp will contain increasing negative values ending at zero. */
167 :     /* Later we add the negative of bp[0] to all values. */
168 :     bp[n] = 0;
169 :     for (j = n - 1; j >= 0; j--) { /* columns in reverse order */
170 :     int i, i1 = ap[j], i2 = ap[j + 1], k, nr;
171 :     if (i1 < i2) AZERO(ind, n);
172 :     for (i = i1; i < i2; i++) {
173 :     ind[ai[i]] = 1;
174 :     for (k = -bp[ai[i] + 1]; k < -bp[ai[i]]; k++) ind[ri[k]] = 1;
175 :     }
176 :     for (k = 0, nr = 0; k < n; k++) if (ind[k]) nr++;
177 :     if ((nr - bp[j + 1]) > bnz) {
178 :     while (nr > (bnz + bp[j + 1])) bnz *= 2;
179 :     ri = Realloc(ri, bnz, int);
180 :     }
181 :     bp[j] = bp[j + 1] - nr;
182 :     for (k = 0, i = -bp[j + 1]; k < n; k++) if (ind[k]) ri[i++] = k;
183 :     }
184 :     bnz = -bp[0];
185 :     bi = INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, bnz));
186 :     bx = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, bnz));
187 :     for (j = 0; j < n; j++) {
188 :     int bpnew = bp[j] + bnz;
189 :     Memcpy(bi + bpnew, ri - bp[j], bp[j + 1] - bp[j]);
190 :     bp[j] = bpnew;
191 :     }
192 :     /* insert code for calculating the actual values here */
193 : bates 1251 for (j = 0; j < bnz; j++) bx[j] = 1;
194 : bates 1248
195 :     Free(ind); Free(ri);
196 :     UNPROTECT(1);
197 :     return ans;
198 :     }
199 : bates 1251 #else
200 :     SEXP dtCMatrix_solve(SEXP a)
201 :     {
202 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
203 :     cs *A = Matrix_as_cs(a);
204 :     int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),
205 :     lo = uplo_P(a)[0] == 'L',
206 :     bnz = 10 * A->n; /* initial estimate of nnz in b */
207 :     int *ti = Calloc(bnz, int), i, j, nz, pos = 0;
208 :     double *tx = Calloc(bnz, double), *wrk = Calloc(A->n, double);
209 :    
210 :     SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
211 :     SET_SLOT(ans, Matrix_DimNamesSym,
212 :     duplicate(GET_SLOT(a, Matrix_DimNamesSym)));
213 :     SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
214 :     SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
215 :     bp[0] = 0;
216 :     for (j = 0; j < A->n; j++) {
217 :     AZERO(wrk, A->n);
218 :     wrk[j] = 1;
219 :     lo ? cs_lsolve(A, wrk) : cs_usolve(A, wrk);
220 :     for (i = 0, nz = 0; i < A->n; i++) if (wrk[i]) nz++;
221 :     bp[j + 1] = nz + bp[j];
222 :     if (bp[j + 1] > bnz) {
223 :     while (bp[j + 1] > bnz) bnz *= 2;
224 :     ti = Realloc(ti, bnz, int);
225 :     tx = Realloc(tx, bnz, double);
226 :     }
227 :     for (i = 0; i < A->n; i++)
228 :     if (wrk[i]) {ti[pos] = i; tx[pos] = wrk[i]; pos++;}
229 :     }
230 :     nz = bp[A->n];
231 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
232 :     Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
233 :    
234 :     Free(A); Free(ti); Free(tx);
235 :     UNPROTECT(1);
236 :     return ans;
237 :     }
238 :     #endif
239 :    
240 :     SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
241 :     {
242 :     int cl = asLogical(classed);
243 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
244 :     cs *A = Matrix_as_cs(a);
245 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
246 :     *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
247 :     getAttrib(b, R_DimSymbol));
248 :     int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
249 :     double *bx;
250 :    
251 :     if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])
252 :     error(_("Dimensions of system to be solved are inconsistent"));
253 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
254 :     /* copy dimnames or Dimnames as well */
255 :     bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),
256 :     REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
257 :     for (j = 0; j < nrhs; j++)
258 :     lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
259 :     Free(A);
260 :     UNPROTECT(1);
261 :     return ans;
262 :     }
263 : bates 1262
264 :     SEXP dtCMatrix_upper_solve(SEXP a)
265 :     {
266 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
267 :     int lo = uplo_P(a)[0] == 'L', unit = diag_P(a)[0] == 'U',
268 :     n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0], nnz,
269 :     *ai = INTEGER(GET_SLOT(a,Matrix_iSym)),
270 :     *ap = INTEGER(GET_SLOT(a, Matrix_pSym)), *bi,
271 :     *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
272 :     int bnz = 10 * ap[n]; /* initial estimate of nnz in b */
273 :     int *ti = Calloc(bnz, int), j, nz;
274 :     double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *tx = Calloc(bnz, double),
275 :     *tmp = Calloc(n, double);
276 :    
277 :     if (lo || (!unit)) error(_("Code written for unit upper triangular unit matrices"));
278 :     bp[0] = 0;
279 :     for (j = 0; j < n; j++) {
280 :     int i, i1 = ap[j + 1];
281 :     AZERO(tmp, n);
282 :     for (i = ap[j]; i < i1; i++) {
283 :     int ii = ai[i], k;
284 :     if (unit) tmp[ii] -= ax[i];
285 :     for (k = bp[ii]; k < bp[ii + 1]; k++) tmp[ti[k]] -= ax[i] * tx[k];
286 :     }
287 :     for (i = 0, nz = 0; i < n; i++) if (tmp[i]) nz++;
288 :     bp[j + 1] = bp[j] + nz;
289 :     if (bp[j + 1] > bnz) {
290 :     while (bp[j + 1] > bnz) bnz *= 2;
291 :     ti = Realloc(ti, bnz, int);
292 :     tx = Realloc(tx, bnz, double);
293 :     }
294 :     i1 = bp[j];
295 :     for (i = 0; i < n; i++) if (tmp[i]) {ti[i1] = i; tx[i1] = tmp[i]; i1++;}
296 :     }
297 :     nz = bp[n];
298 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
299 :     Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
300 :     Free(tmp); Free(tx); Free(ti);
301 :     SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
302 :     SET_SLOT(ans, Matrix_DimNamesSym, duplicate(GET_SLOT(a, Matrix_DimNamesSym)));
303 :     SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
304 :     SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
305 :     UNPROTECT(1);
306 :     return ans;
307 :     }

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