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 1248 - (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 10
4 :     SEXP tsc_validate(SEXP x)
5 :     {
6 : maechler 890 return triangularMatrix_validate(x);
7 : maechler 895 /* see ./dsCMatrix.c or ./dtpMatrix.c on how to do more testing here */
8 : bates 10 }
9 :    
10 : bates 1209 #if 0
11 : bates 10 SEXP tsc_transpose(SEXP x)
12 :     {
13 : bates 1209 cholmod_sparse *cx = as_cholmod_sparse(x);
14 :    
15 : maechler 945 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix"))),
16 : bates 10 islot = GET_SLOT(x, Matrix_iSym);
17 :     int nnz = length(islot),
18 : bates 338 *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
19 : maechler 951 int up = uplo_P(x)[0] == 'U';
20 : bates 10
21 : bates 588 adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
22 : bates 10 adims[0] = xdims[1]; adims[1] = xdims[0];
23 : maechler 945
24 : maechler 951 if(*diag_P(x) == 'U')
25 : maechler 945 SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(x, Matrix_diagSym)));
26 : bates 741 SET_SLOT(ans, Matrix_uploSym, mkString(up ? "L" : "U"));
27 : maechler 945
28 : bates 588 csc_compTr(xdims[0], xdims[1], nnz,
29 :     INTEGER(GET_SLOT(x, Matrix_pSym)), INTEGER(islot),
30 :     REAL(GET_SLOT(x, Matrix_xSym)),
31 :     INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, xdims[0] + 1)),
32 :     INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)),
33 :     REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)));
34 : bates 10 UNPROTECT(1);
35 :     return ans;
36 :     }
37 : bates 1209 #endif
38 : bates 70
39 : bates 478 SEXP tsc_to_dgTMatrix(SEXP x)
40 : bates 70 {
41 :     SEXP ans;
42 : maechler 951 if (*diag_P(x) != 'U')
43 : bates 677 ans = compressed_to_dgTMatrix(x, ScalarLogical(1));
44 : bates 70 else { /* unit triangular matrix */
45 : maechler 534 SEXP islot = GET_SLOT(x, Matrix_iSym),
46 : bates 70 pslot = GET_SLOT(x, Matrix_pSym);
47 :     int *ai, *aj, j,
48 :     n = length(pslot) - 1,
49 :     nod = length(islot),
50 :     nout = n + nod,
51 :     *p = INTEGER(pslot);
52 :     double *ax;
53 : maechler 534
54 : bates 478 ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix")));
55 : bates 70 SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
56 :     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nout));
57 :     ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
58 :     Memcpy(ai, INTEGER(islot), nod);
59 :     SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, nout));
60 :     aj = INTEGER(GET_SLOT(ans, Matrix_jSym));
61 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nout));
62 :     ax = REAL(GET_SLOT(ans, Matrix_xSym));
63 :     Memcpy(ax, REAL(GET_SLOT(x, Matrix_xSym)), nod);
64 :     for (j = 0; j < n; j++) {
65 :     int jj, npj = nod + j, p2 = p[j+1];
66 :     aj[npj] = ai[npj] = j;
67 :     ax[npj] = 1.;
68 :     for (jj = p[j]; jj < p2; jj++) aj[jj] = j;
69 :     }
70 :     UNPROTECT(1);
71 :     }
72 :     return ans;
73 :     }
74 : bates 338
75 : maechler 534 /**
76 : bates 358 * Derive the column pointer vector for the inverse of L from the parent array
77 : maechler 534 *
78 : bates 358 * @param n length of parent array
79 :     * @param countDiag 0 for a unit triangular matrix with implicit diagonal, otherwise 1
80 :     * @param pr parent vector describing the elimination tree
81 :     * @param ap array of length n+1 to be filled with the column pointers
82 : maechler 534 *
83 : bates 358 * @return the number of non-zero entries (ap[n])
84 :     */
85 :     int parent_inv_ap(int n, int countDiag, const int pr[], int ap[])
86 :     {
87 :     int *sz = Calloc(n, int), j;
88 :    
89 :     for (j = n - 1; j >= 0; j--) {
90 :     int parent = pr[j];
91 :     sz[j] = (parent < 0) ? countDiag : (1 + sz[parent]);
92 :     }
93 :     ap[0] = 0;
94 :     for (j = 0; j < n; j++)
95 :     ap[j+1] = ap[j] + sz[j];
96 :     Free(sz);
97 :     return ap[n];
98 :     }
99 :    
100 : maechler 534 /**
101 : bates 358 * Derive the row index array for the inverse of L from the parent array
102 : maechler 534 *
103 : bates 358 * @param n length of parent array
104 :     * @param countDiag 0 for a unit triangular matrix with implicit diagonal, otherwise 1
105 :     * @param pr parent vector describing the elimination tree
106 :     * @param ai row index vector of length ap[n]
107 :     */
108 :     void parent_inv_ai(int n, int countDiag, const int pr[], int ai[])
109 :     {
110 :     int j, k, pos = 0;
111 :     for (j = 0; j < n; j++) {
112 :     if (countDiag) ai[pos++] = j;
113 :     for (k = pr[j]; k >= 0; k = pr[k]) ai[pos++] = k;
114 :     }
115 :     }
116 : maechler 534
117 : bates 338 SEXP Parent_inverse(SEXP par, SEXP unitdiag)
118 :     {
119 : bates 478 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
120 : bates 338 int *ap, *ai, *dims, *pr = INTEGER(par),
121 :     countDiag = 1 - asLogical(unitdiag),
122 : bates 367 j, n = length(par), nnz;
123 : bates 338 double *ax;
124 : maechler 534
125 : bates 582 if (!isInteger(par)) error(_("par argument must be an integer vector"));
126 : bates 338 SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, n + 1));
127 :     ap = INTEGER(GET_SLOT(ans, Matrix_pSym));
128 : bates 358 nnz = parent_inv_ap(n, countDiag, pr, ap);
129 : bates 338 SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
130 :     ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
131 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
132 :     ax = REAL(GET_SLOT(ans, Matrix_xSym));
133 :     for (j = 0; j < nnz; j++) ax[j] = 1.;
134 :     SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
135 :     dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
136 :     dims[0] = dims[1] = n;
137 : maechler 534 SET_SLOT(ans, Matrix_uploSym, mkString("L"));
138 :     SET_SLOT(ans, Matrix_diagSym, (countDiag ? mkString("N") : mkString("U")));
139 : bates 358 parent_inv_ai(n, countDiag, pr, ai);
140 : bates 338 UNPROTECT(1);
141 :     return ans;
142 :     }
143 : bates 1248
144 :     SEXP dtCMatrix_solve(SEXP a)
145 :     {
146 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
147 :     int lo = uplo_P(a)[0] == 'L', unit = diag_P(a)[0] == 'U',
148 :     n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0], nnz,
149 :     *ai = INTEGER(GET_SLOT(a,Matrix_iSym)),
150 :     *ap = INTEGER(GET_SLOT(a, Matrix_pSym)), *bi,
151 :     *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
152 :     int bnz = 10 * ap[n]; /* initial estimate of nnz in b */
153 :     int *ri = Calloc(bnz, int), *ind = Calloc(n, int), j;
154 :     double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *bx;
155 :    
156 :     SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
157 :     SET_SLOT(ans, Matrix_DimNamesSym,
158 :     duplicate(GET_SLOT(a, Matrix_DimNamesSym)));
159 :     SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
160 :     SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
161 :    
162 :     if (!(lo && unit))
163 :     error("code for non-unit or upper triangular not yet written");
164 :     /* Initially bp will contain increasing negative values ending at zero. */
165 :     /* Later we add the negative of bp[0] to all values. */
166 :     bp[n] = 0;
167 :     for (j = n - 1; j >= 0; j--) { /* columns in reverse order */
168 :     int i, i1 = ap[j], i2 = ap[j + 1], k, nr;
169 :     if (i1 < i2) AZERO(ind, n);
170 :     for (i = i1; i < i2; i++) {
171 :     ind[ai[i]] = 1;
172 :     for (k = -bp[ai[i] + 1]; k < -bp[ai[i]]; k++) ind[ri[k]] = 1;
173 :     }
174 :     for (k = 0, nr = 0; k < n; k++) if (ind[k]) nr++;
175 :     if ((nr - bp[j + 1]) > bnz) {
176 :     while (nr > (bnz + bp[j + 1])) bnz *= 2;
177 :     ri = Realloc(ri, bnz, int);
178 :     }
179 :     bp[j] = bp[j + 1] - nr;
180 :     for (k = 0, i = -bp[j + 1]; k < n; k++) if (ind[k]) ri[i++] = k;
181 :     }
182 :     bnz = -bp[0];
183 :     bi = INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, bnz));
184 :     bx = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, bnz));
185 :     for (j = 0; j < n; j++) {
186 :     int bpnew = bp[j] + bnz;
187 :     Memcpy(bi + bpnew, ri - bp[j], bp[j + 1] - bp[j]);
188 :     bp[j] = bpnew;
189 :     }
190 :     /* insert code for calculating the actual values here */
191 :     for (j = 0; j < bnz; j++) bx[j] = NA_REAL;
192 :    
193 :     Free(ind); Free(ri);
194 :     UNPROTECT(1);
195 :     return ans;
196 :     }

root@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