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 1968 - (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 : maechler 1968 #define RETURN(_CH_) UNPROTECT(1); return (_CH_);
6 :    
7 :     /* This is used for *BOTH* triangular and symmetric Csparse: */
8 : maechler 1710 SEXP tCMatrix_validate(SEXP x)
9 :     {
10 :     SEXP val = xCMatrix_validate(x);/* checks x slot */
11 :     if(isString(val))
12 :     return(val);
13 :     else {
14 :     SEXP
15 :     islot = GET_SLOT(x, Matrix_iSym),
16 :     pslot = GET_SLOT(x, Matrix_pSym);
17 :     int uploT = (*uplo_P(x) == 'U'),
18 :     k, nnz = length(islot),
19 :     *xi = INTEGER(islot),
20 :     *xj = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
21 :    
22 :     expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xj);
23 :    
24 :     /* Maybe FIXME: ">" should be ">=" for diag = 'U' (uplo = 'U') */
25 :     if(uploT) {
26 :     for (k = 0; k < nnz; k++)
27 :     if(xi[k] > xj[k]) {
28 :     RETURN(mkString(_("uplo='U' must not have sparse entries in lower diagonal")));
29 :     }
30 :     }
31 :     else {
32 :     for (k = 0; k < nnz; k++)
33 :     if(xi[k] < xj[k]) {
34 :     RETURN(mkString(_("uplo='L' must not have sparse entries in upper diagonal")));
35 :     }
36 :     }
37 :    
38 :     RETURN(ScalarLogical(1));
39 :     }
40 :     }
41 : maechler 1968
42 :     /* This is used for *BOTH* triangular and symmetric Rsparse: */
43 :     SEXP tRMatrix_validate(SEXP x)
44 :     {
45 :     SEXP val = xRMatrix_validate(x);/* checks x slot */
46 :     if(isString(val))
47 :     return(val);
48 :     else {
49 :     SEXP
50 :     jslot = GET_SLOT(x, Matrix_jSym),
51 :     pslot = GET_SLOT(x, Matrix_pSym);
52 :     int uploT = (*uplo_P(x) == 'U'),
53 :     k, nnz = length(jslot),
54 :     *xj = INTEGER(jslot),
55 :     *xi = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
56 :    
57 :     expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xi);
58 :    
59 :     /* Maybe FIXME: ">" should be ">=" for diag = 'U' (uplo = 'U') */
60 :     if(uploT) {
61 :     for (k = 0; k < nnz; k++)
62 :     if(xi[k] > xj[k]) {
63 :     RETURN(mkString(_("uplo='U' must not have sparse entries in lower diagonal")));
64 :     }
65 :     }
66 :     else {
67 :     for (k = 0; k < nnz; k++)
68 :     if(xi[k] < xj[k]) {
69 :     RETURN(mkString(_("uplo='L' must not have sparse entries in upper diagonal")));
70 :     }
71 :     }
72 :    
73 :     RETURN(ScalarLogical(1));
74 :     }
75 :     }
76 :    
77 :    
78 : maechler 1710 #undef RETURN
79 :    
80 : maechler 534 /**
81 : bates 358 * Derive the column pointer vector for the inverse of L from the parent array
82 : maechler 534 *
83 : bates 358 * @param n length of parent array
84 :     * @param countDiag 0 for a unit triangular matrix with implicit diagonal, otherwise 1
85 :     * @param pr parent vector describing the elimination tree
86 :     * @param ap array of length n+1 to be filled with the column pointers
87 : maechler 534 *
88 : bates 358 * @return the number of non-zero entries (ap[n])
89 :     */
90 :     int parent_inv_ap(int n, int countDiag, const int pr[], int ap[])
91 :     {
92 : maechler 1960 int *sz = Alloca(n, int), j;
93 :     R_CheckStack();
94 : bates 358
95 :     for (j = n - 1; j >= 0; j--) {
96 :     int parent = pr[j];
97 :     sz[j] = (parent < 0) ? countDiag : (1 + sz[parent]);
98 :     }
99 :     ap[0] = 0;
100 :     for (j = 0; j < n; j++)
101 :     ap[j+1] = ap[j] + sz[j];
102 :     return ap[n];
103 :     }
104 :    
105 : maechler 534 /**
106 : bates 358 * Derive the row index array for the inverse of L from the parent array
107 : maechler 534 *
108 : bates 358 * @param n length of parent array
109 :     * @param countDiag 0 for a unit triangular matrix with implicit diagonal, otherwise 1
110 :     * @param pr parent vector describing the elimination tree
111 :     * @param ai row index vector of length ap[n]
112 :     */
113 :     void parent_inv_ai(int n, int countDiag, const int pr[], int ai[])
114 :     {
115 :     int j, k, pos = 0;
116 :     for (j = 0; j < n; j++) {
117 :     if (countDiag) ai[pos++] = j;
118 :     for (k = pr[j]; k >= 0; k = pr[k]) ai[pos++] = k;
119 :     }
120 :     }
121 : maechler 534
122 : bates 338 SEXP Parent_inverse(SEXP par, SEXP unitdiag)
123 :     {
124 : bates 478 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
125 : bates 338 int *ap, *ai, *dims, *pr = INTEGER(par),
126 :     countDiag = 1 - asLogical(unitdiag),
127 : bates 367 j, n = length(par), nnz;
128 : bates 338 double *ax;
129 : maechler 534
130 : bates 582 if (!isInteger(par)) error(_("par argument must be an integer vector"));
131 : bates 338 SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, n + 1));
132 :     ap = INTEGER(GET_SLOT(ans, Matrix_pSym));
133 : bates 358 nnz = parent_inv_ap(n, countDiag, pr, ap);
134 : bates 338 SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
135 :     ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
136 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
137 :     ax = REAL(GET_SLOT(ans, Matrix_xSym));
138 :     for (j = 0; j < nnz; j++) ax[j] = 1.;
139 :     SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
140 :     dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
141 :     dims[0] = dims[1] = n;
142 : maechler 534 SET_SLOT(ans, Matrix_uploSym, mkString("L"));
143 :     SET_SLOT(ans, Matrix_diagSym, (countDiag ? mkString("N") : mkString("U")));
144 : bates 358 parent_inv_ai(n, countDiag, pr, ai);
145 : bates 338 UNPROTECT(1);
146 :     return ans;
147 :     }
148 : bates 1248
149 :     SEXP dtCMatrix_solve(SEXP a)
150 :     {
151 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
152 : maechler 1960 CSP A = AS_CSP(a);
153 : bates 1251 int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),
154 : maechler 1960 bnz = 10 * A->n, /* initial estimate of nnz in b */
155 :     lo = uplo_P(a)[0] == 'L', top;
156 :     /* These arrays must use Calloc because of possible Realloc */
157 : maechler 1798 int *ti = Calloc(bnz, int), p, j, nz, pos = 0;
158 : maechler 1960 double *tx = Calloc(bnz, double);
159 : maechler 1798 cs *u = cs_spalloc(A->n, 1,1,1,0); /* Sparse unit vector */
160 : maechler 1960 double *wrk = Alloca(A->n, double);
161 :     int *xi = Alloca(2*A->n, int); /* for cs_reach */
162 :     R_CheckStack();
163 : bates 1251
164 : maechler 1968 slot_dup(ans, a, Matrix_DimSym);
165 : maechler 1736 SET_DimNames(ans, a);
166 : maechler 1968 slot_dup(ans, a, Matrix_uploSym);
167 :     slot_dup(ans, a, Matrix_diagSym);
168 : maechler 1798 /* initialize the "constant part" of the sparse unit vector */
169 :     u->x[0] = 1.;
170 :     u->p[0] = 0; u->p[1] = 1;
171 : bates 1251 bp[0] = 0;
172 :     for (j = 0; j < A->n; j++) {
173 : maechler 1798 u->i[0] = j; /* u := j'th unit vector */
174 :     /* (wrk[top:n],xi[top:n]) := A^{-1} u : */
175 :     top = cs_spsolve (A, u, 0, xi, wrk, 0, lo);
176 :     nz = A->n - top;
177 : bates 1251 bp[j + 1] = nz + bp[j];
178 :     if (bp[j + 1] > bnz) {
179 :     while (bp[j + 1] > bnz) bnz *= 2;
180 :     ti = Realloc(ti, bnz, int);
181 :     tx = Realloc(tx, bnz, double);
182 :     }
183 : maechler 1798 if (lo)
184 :     for(p = top; p < A->n; p++, pos++) {
185 :     ti[pos] = xi[p];
186 :     tx[pos] = wrk[xi[p]];
187 :     }
188 :     else /* upper triagonal */
189 :     for(p = A->n - 1; p >= top; p--, pos++) {
190 :     ti[pos] = xi[p];
191 :     tx[pos] = wrk[xi[p]];
192 :     }
193 : bates 1251 }
194 :     nz = bp[A->n];
195 : maechler 1798 Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
196 :     Memcpy( REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
197 : bates 1251
198 : maechler 1960 Free(ti); Free(tx); cs_spfree(u);
199 : bates 1251 UNPROTECT(1);
200 :     return ans;
201 :     }
202 :    
203 :     SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
204 :     {
205 :     int cl = asLogical(classed);
206 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
207 : maechler 1960 CSP A = AS_CSP(a);
208 : bates 1251 int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
209 :     *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
210 :     getAttrib(b, R_DimSymbol));
211 :     int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
212 :     double *bx;
213 : maechler 1960 R_CheckStack();
214 : bates 1251
215 :     if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])
216 :     error(_("Dimensions of system to be solved are inconsistent"));
217 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
218 : maechler 1736 /* FIXME: copy dimnames or Dimnames as well */
219 : bates 1251 bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),
220 :     REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
221 :     for (j = 0; j < nrhs; j++)
222 :     lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
223 :     UNPROTECT(1);
224 :     return ans;
225 :     }
226 : bates 1262
227 :     SEXP dtCMatrix_upper_solve(SEXP a)
228 :     {
229 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
230 :     int lo = uplo_P(a)[0] == 'L', unit = diag_P(a)[0] == 'U',
231 : bates 1310 n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0],
232 : bates 1262 *ai = INTEGER(GET_SLOT(a,Matrix_iSym)),
233 : maechler 1661 *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),
234 : bates 1262 *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
235 :     int bnz = 10 * ap[n]; /* initial estimate of nnz in b */
236 : maechler 1960 int *ti = Alloca(bnz, int), j, nz;
237 :     double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *tx = Alloca(bnz, double),
238 :     *tmp = Alloca(n, double);
239 :     R_CheckStack();
240 : maechler 1661
241 : bates 1386 if (lo || (!unit))
242 :     error(_("Code written for unit upper triangular unit matrices"));
243 : bates 1262 bp[0] = 0;
244 :     for (j = 0; j < n; j++) {
245 :     int i, i1 = ap[j + 1];
246 :     AZERO(tmp, n);
247 :     for (i = ap[j]; i < i1; i++) {
248 :     int ii = ai[i], k;
249 :     if (unit) tmp[ii] -= ax[i];
250 :     for (k = bp[ii]; k < bp[ii + 1]; k++) tmp[ti[k]] -= ax[i] * tx[k];
251 :     }
252 :     for (i = 0, nz = 0; i < n; i++) if (tmp[i]) nz++;
253 :     bp[j + 1] = bp[j] + nz;
254 :     if (bp[j + 1] > bnz) {
255 :     while (bp[j + 1] > bnz) bnz *= 2;
256 :     ti = Realloc(ti, bnz, int);
257 :     tx = Realloc(tx, bnz, double);
258 :     }
259 :     i1 = bp[j];
260 :     for (i = 0; i < n; i++) if (tmp[i]) {ti[i1] = i; tx[i1] = tmp[i]; i1++;}
261 :     }
262 :     nz = bp[n];
263 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
264 :     Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
265 : maechler 1968 slot_dup(ans, a, Matrix_DimSym);
266 : maechler 1736 SET_DimNames(ans, a);
267 : maechler 1968 slot_dup(ans, a, Matrix_uploSym);
268 :     slot_dup(ans, a, Matrix_diagSym);
269 : bates 1262 UNPROTECT(1);
270 :     return ans;
271 :     }

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