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

R-Forge@R-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business University of Wisconsin - Madison Powered By FusionForge