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 1507 - (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 : bates 1507 SEXP dtCMatrix_validate(SEXP x)
6 : bates 10 {
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 1507 #if 0 /* no longer used */
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 1507 #endif
48 : bates 338
49 : maechler 534 /**
50 : bates 358 * Derive the column pointer vector for the inverse of L from the parent array
51 : maechler 534 *
52 : bates 358 * @param n length of parent array
53 :     * @param countDiag 0 for a unit triangular matrix with implicit diagonal, otherwise 1
54 :     * @param pr parent vector describing the elimination tree
55 :     * @param ap array of length n+1 to be filled with the column pointers
56 : maechler 534 *
57 : bates 358 * @return the number of non-zero entries (ap[n])
58 :     */
59 :     int parent_inv_ap(int n, int countDiag, const int pr[], int ap[])
60 :     {
61 :     int *sz = Calloc(n, int), j;
62 :    
63 :     for (j = n - 1; j >= 0; j--) {
64 :     int parent = pr[j];
65 :     sz[j] = (parent < 0) ? countDiag : (1 + sz[parent]);
66 :     }
67 :     ap[0] = 0;
68 :     for (j = 0; j < n; j++)
69 :     ap[j+1] = ap[j] + sz[j];
70 :     Free(sz);
71 :     return ap[n];
72 :     }
73 :    
74 : maechler 534 /**
75 : bates 358 * Derive the row index array for the inverse of L from the parent array
76 : maechler 534 *
77 : bates 358 * @param n length of parent array
78 :     * @param countDiag 0 for a unit triangular matrix with implicit diagonal, otherwise 1
79 :     * @param pr parent vector describing the elimination tree
80 :     * @param ai row index vector of length ap[n]
81 :     */
82 :     void parent_inv_ai(int n, int countDiag, const int pr[], int ai[])
83 :     {
84 :     int j, k, pos = 0;
85 :     for (j = 0; j < n; j++) {
86 :     if (countDiag) ai[pos++] = j;
87 :     for (k = pr[j]; k >= 0; k = pr[k]) ai[pos++] = k;
88 :     }
89 :     }
90 : maechler 534
91 : bates 338 SEXP Parent_inverse(SEXP par, SEXP unitdiag)
92 :     {
93 : bates 478 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
94 : bates 338 int *ap, *ai, *dims, *pr = INTEGER(par),
95 :     countDiag = 1 - asLogical(unitdiag),
96 : bates 367 j, n = length(par), nnz;
97 : bates 338 double *ax;
98 : maechler 534
99 : bates 582 if (!isInteger(par)) error(_("par argument must be an integer vector"));
100 : bates 338 SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, n + 1));
101 :     ap = INTEGER(GET_SLOT(ans, Matrix_pSym));
102 : bates 358 nnz = parent_inv_ap(n, countDiag, pr, ap);
103 : bates 338 SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
104 :     ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
105 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
106 :     ax = REAL(GET_SLOT(ans, Matrix_xSym));
107 :     for (j = 0; j < nnz; j++) ax[j] = 1.;
108 :     SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
109 :     dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
110 :     dims[0] = dims[1] = n;
111 : maechler 534 SET_SLOT(ans, Matrix_uploSym, mkString("L"));
112 :     SET_SLOT(ans, Matrix_diagSym, (countDiag ? mkString("N") : mkString("U")));
113 : bates 358 parent_inv_ai(n, countDiag, pr, ai);
114 : bates 338 UNPROTECT(1);
115 :     return ans;
116 :     }
117 : bates 1248
118 :     SEXP dtCMatrix_solve(SEXP a)
119 :     {
120 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
121 : bates 1251 cs *A = Matrix_as_cs(a);
122 :     int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),
123 :     lo = uplo_P(a)[0] == 'L',
124 :     bnz = 10 * A->n; /* initial estimate of nnz in b */
125 :     int *ti = Calloc(bnz, int), i, j, nz, pos = 0;
126 :     double *tx = Calloc(bnz, double), *wrk = Calloc(A->n, double);
127 :    
128 :     SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
129 :     SET_SLOT(ans, Matrix_DimNamesSym,
130 :     duplicate(GET_SLOT(a, Matrix_DimNamesSym)));
131 :     SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
132 :     SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
133 :     bp[0] = 0;
134 :     for (j = 0; j < A->n; j++) {
135 :     AZERO(wrk, A->n);
136 :     wrk[j] = 1;
137 :     lo ? cs_lsolve(A, wrk) : cs_usolve(A, wrk);
138 :     for (i = 0, nz = 0; i < A->n; i++) if (wrk[i]) nz++;
139 :     bp[j + 1] = nz + bp[j];
140 :     if (bp[j + 1] > bnz) {
141 :     while (bp[j + 1] > bnz) bnz *= 2;
142 :     ti = Realloc(ti, bnz, int);
143 :     tx = Realloc(tx, bnz, double);
144 :     }
145 :     for (i = 0; i < A->n; i++)
146 :     if (wrk[i]) {ti[pos] = i; tx[pos] = wrk[i]; pos++;}
147 :     }
148 :     nz = bp[A->n];
149 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
150 :     Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
151 :    
152 :     Free(A); Free(ti); Free(tx);
153 :     UNPROTECT(1);
154 :     return ans;
155 :     }
156 :    
157 :     SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
158 :     {
159 :     int cl = asLogical(classed);
160 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
161 :     cs *A = Matrix_as_cs(a);
162 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
163 :     *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
164 :     getAttrib(b, R_DimSymbol));
165 :     int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
166 :     double *bx;
167 :    
168 :     if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])
169 :     error(_("Dimensions of system to be solved are inconsistent"));
170 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
171 :     /* copy dimnames or Dimnames as well */
172 :     bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),
173 :     REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
174 :     for (j = 0; j < nrhs; j++)
175 :     lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
176 :     Free(A);
177 :     UNPROTECT(1);
178 :     return ans;
179 :     }
180 : bates 1262
181 :     SEXP dtCMatrix_upper_solve(SEXP a)
182 :     {
183 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
184 :     int lo = uplo_P(a)[0] == 'L', unit = diag_P(a)[0] == 'U',
185 : bates 1310 n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0],
186 : bates 1262 *ai = INTEGER(GET_SLOT(a,Matrix_iSym)),
187 : bates 1310 *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),
188 : bates 1262 *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
189 :     int bnz = 10 * ap[n]; /* initial estimate of nnz in b */
190 :     int *ti = Calloc(bnz, int), j, nz;
191 :     double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *tx = Calloc(bnz, double),
192 :     *tmp = Calloc(n, double);
193 :    
194 : bates 1386 if (lo || (!unit))
195 :     error(_("Code written for unit upper triangular unit matrices"));
196 : bates 1262 bp[0] = 0;
197 :     for (j = 0; j < n; j++) {
198 :     int i, i1 = ap[j + 1];
199 :     AZERO(tmp, n);
200 :     for (i = ap[j]; i < i1; i++) {
201 :     int ii = ai[i], k;
202 :     if (unit) tmp[ii] -= ax[i];
203 :     for (k = bp[ii]; k < bp[ii + 1]; k++) tmp[ti[k]] -= ax[i] * tx[k];
204 :     }
205 :     for (i = 0, nz = 0; i < n; i++) if (tmp[i]) nz++;
206 :     bp[j + 1] = bp[j] + nz;
207 :     if (bp[j + 1] > bnz) {
208 :     while (bp[j + 1] > bnz) bnz *= 2;
209 :     ti = Realloc(ti, bnz, int);
210 :     tx = Realloc(tx, bnz, double);
211 :     }
212 :     i1 = bp[j];
213 :     for (i = 0; i < n; i++) if (tmp[i]) {ti[i1] = i; tx[i1] = tmp[i]; i1++;}
214 :     }
215 :     nz = bp[n];
216 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
217 :     Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
218 :     Free(tmp); Free(tx); Free(ti);
219 :     SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
220 :     SET_SLOT(ans, Matrix_DimNamesSym, duplicate(GET_SLOT(a, Matrix_DimNamesSym)));
221 :     SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
222 :     SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
223 :     UNPROTECT(1);
224 :     return ans;
225 :     }

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