SCM

SCM Repository

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

Annotation of /pkg/src/dsCMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 683 - (view) (download) (as text)

1 : bates 478 #include "dsCMatrix.h"
2 : bates 10
3 : bates 478 SEXP dsCMatrix_validate(SEXP obj)
4 : bates 10 {
5 : bates 592 SEXP val = check_scalar_string(GET_SLOT(obj, Matrix_uploSym),
6 :     "LU", "uplo");
7 : bates 10 int *Dim = INTEGER(GET_SLOT(obj, Matrix_DimSym));
8 : maechler 534
9 : bates 592 if (isString(val)) return val;
10 : bates 10 if (Dim[0] != Dim[1])
11 : bates 582 return mkString(_("Symmetric matrix must be square"));
12 : bates 10 csc_check_column_sorting(obj);
13 :     return ScalarLogical(1);
14 :     }
15 :    
16 : bates 478 SEXP dsCMatrix_chol(SEXP x, SEXP pivot)
17 : bates 10 {
18 : bates 257 SEXP pSlot = GET_SLOT(x, Matrix_pSym), xorig = x;
19 :     int *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
20 :     *Ap = INTEGER(pSlot),
21 :     *Lp, *Parent, info,
22 : maechler 534 lo = CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'L',
23 : bates 257 n = length(pSlot)-1,
24 :     nnz, piv = asLogical(pivot);
25 : bates 488 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dCholCMatrix")));
26 : bates 492 int *P, *Pinv;
27 : bates 257 double *Ax;
28 : bates 10
29 : bates 492 /* FIXME: Check if there is a Cholesky factorization. If yes,
30 :     check if the permutation status matches that of the call. If
31 :     so, return it. */
32 : maechler 534
33 : bates 257 if (lo) {
34 :     x = PROTECT(ssc_transpose(x));
35 :     Ai = INTEGER(GET_SLOT(x, Matrix_iSym));
36 :     Ap = INTEGER(GET_SLOT(x, Matrix_pSym));
37 :     }
38 : bates 492 SET_SLOT(val, Matrix_uploSym, mkString("L"));
39 :     SET_SLOT(val, Matrix_diagSym, mkString("U"));
40 : bates 257 SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
41 :     SET_SLOT(val, Matrix_ParentSym, allocVector(INTSXP, n));
42 :     Parent = INTEGER(GET_SLOT(val, Matrix_ParentSym));
43 :     SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, n + 1));
44 :     Lp = INTEGER(GET_SLOT(val, Matrix_pSym));
45 : bates 492 SET_SLOT(val, Matrix_permSym, allocVector(INTSXP, n));
46 :     P = INTEGER(GET_SLOT(val, Matrix_permSym));
47 : bates 10 if (piv) {
48 : bates 478 SEXP trip = PROTECT(dsCMatrix_to_dgTMatrix(x));
49 : bates 257 SEXP Ti = GET_SLOT(trip, Matrix_iSym);
50 : bates 10
51 : bates 257 /* determine the permutation with Metis */
52 :     Pinv = Calloc(n, int);
53 :     ssc_metis_order(n, Ap, Ai, P, Pinv);
54 :     /* create a symmetrized form of x */
55 :     nnz = length(Ti);
56 :     Ai = Calloc(nnz, int);
57 :     Ax = Calloc(nnz, double);
58 :     Ap = Calloc(n + 1, int);
59 : bates 488 triplet_to_col(n, n, nnz, INTEGER(Ti),
60 : bates 257 INTEGER(GET_SLOT(trip, Matrix_jSym)),
61 :     REAL(GET_SLOT(trip, Matrix_xSym)),
62 :     Ap, Ai, Ax);
63 : bates 492 UNPROTECT(1);
64 :     } else {
65 :     int i;
66 :     for (i = 0; i < n; i++) P[i] = i;
67 :     Ax = REAL(GET_SLOT(x, Matrix_xSym));
68 :    
69 :     }
70 :     R_ldl_symbolic(n, Ap, Ai, Lp, Parent, (piv) ? P : (int *)NULL,
71 :     (piv) ? Pinv : (int *)NULL);
72 : bates 257 nnz = Lp[n];
73 : bates 10 SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));
74 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));
75 : maechler 534 SET_SLOT(val, Matrix_DSym, allocVector(REALSXP, n));
76 : bates 492 info = R_ldl_numeric(n, Ap, Ai, Ax, Lp, Parent,
77 :     INTEGER(GET_SLOT(val, Matrix_iSym)),
78 :     REAL(GET_SLOT(val, Matrix_xSym)),
79 :     REAL(GET_SLOT(val, Matrix_DSym)),
80 :     (piv) ? P : (int *)NULL,
81 :     (piv) ? Pinv : (int *)NULL);
82 : bates 257 if (info != n)
83 : bates 582 error(_("Leading minor of size %d (possibly after permutation) is indefinite"),
84 : bates 257 info + 1);
85 :     if (piv) {
86 :     Free(Pinv); Free(Ax); Free(Ai); Free(Ap);
87 :     }
88 :     UNPROTECT(lo ? 2 : 1);
89 : bates 476 return set_factors(xorig, val, "Cholesky");
90 : bates 10 }
91 :    
92 : bates 683 SEXP dsCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
93 : bates 10 {
94 : bates 683 int cl = asLogical(classed);
95 : bates 476 SEXP Chol = get_factors(a, "Cholesky"), perm,
96 : bates 683 bdP = cl ? GET_SLOT(b, Matrix_DimSym) : getAttrib(b, R_DimSymbol),
97 :     val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
98 : bates 10 int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
99 : bates 683 *bdims = INTEGER(bdP), *Li, *Lp, j, piv;
100 :     int n = adims[1], ncol = bdims[1];
101 :     double *Lx, *D, *in = REAL(cl ? GET_SLOT(b, Matrix_xSym) : b),
102 :     *out = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n * ncol)),
103 :     *tmp = (double *) NULL;
104 : bates 10
105 : bates 683 if (!cl && !(isReal(b) && isMatrix(b)))
106 : bates 582 error(_("Argument b must be a numeric matrix"));
107 : bates 683 if (*adims != *bdims || ncol < 1 || *adims < 1)
108 : bates 582 error(_("Dimensions of system to be solved are inconsistent"));
109 : bates 478 if (Chol == R_NilValue) Chol = dsCMatrix_chol(a, ScalarLogical(1));
110 : bates 257 perm = GET_SLOT(Chol, Matrix_permSym);
111 :     piv = length(perm);
112 :     if (piv) tmp = Calloc(n, double);
113 :     Li = INTEGER(GET_SLOT(Chol, Matrix_iSym));
114 :     Lp = INTEGER(GET_SLOT(Chol, Matrix_pSym));
115 :     Lx = REAL(GET_SLOT(Chol, Matrix_xSym));
116 :     D = REAL(GET_SLOT(Chol, Matrix_DSym));
117 : bates 308 for (j = 0; j < ncol; j++, in += n, out += n) {
118 : bates 366 if (piv) R_ldl_perm(n, out, in, INTEGER(perm));
119 : bates 257 else Memcpy(out, in, n);
120 : bates 366 R_ldl_lsolve(n, out, Lp, Li, Lx);
121 :     R_ldl_dsolve(n, out, D);
122 :     R_ldl_ltsolve(n, out, Lp, Li, Lx);
123 :     if (piv) R_ldl_permt(n, out, Memcpy(tmp, out, n), INTEGER(perm));
124 : bates 10 }
125 : bates 257 if (piv) Free(tmp);
126 : bates 10 UNPROTECT(1);
127 :     return val;
128 :     }
129 :    
130 : bates 478 SEXP dsCMatrix_inverse_factor(SEXP A)
131 : bates 10 {
132 : bates 257 return R_NilValue; /* FIXME: Write this function. */
133 : bates 10 }
134 :    
135 :     SEXP ssc_transpose(SEXP x)
136 :     {
137 : bates 587 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsCMatrix"))),
138 : bates 10 islot = GET_SLOT(x, Matrix_iSym);
139 : bates 587 int nnz = length(islot), *adims,
140 : bates 10 *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
141 :    
142 : bates 587 adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
143 : bates 10 adims[0] = xdims[1]; adims[1] = xdims[0];
144 : maechler 534 if (CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'U')
145 :     SET_SLOT(ans, Matrix_uploSym, mkString("L"));
146 : maechler 614 else
147 :     SET_SLOT(ans, Matrix_uploSym, mkString("U"));
148 : bates 587 csc_compTr(xdims[0], xdims[1], nnz,
149 :     INTEGER(GET_SLOT(x, Matrix_pSym)), INTEGER(islot),
150 :     REAL(GET_SLOT(x, Matrix_xSym)),
151 :     INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, xdims[0] + 1)),
152 :     INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)),
153 :     REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)));
154 : bates 10 UNPROTECT(1);
155 :     return ans;
156 :     }
157 : bates 70
158 : bates 478 SEXP dsCMatrix_to_dgTMatrix(SEXP x)
159 : bates 70 {
160 :     SEXP
161 : bates 478 ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix"))),
162 : bates 70 islot = GET_SLOT(x, Matrix_iSym),
163 :     pslot = GET_SLOT(x, Matrix_pSym);
164 :     int *ai, *aj, *iv = INTEGER(islot),
165 :     j, jj, nnz = length(islot), nout,
166 :     n = length(pslot) - 1,
167 :     *p = INTEGER(pslot), pos;
168 :     double *ax, *xv = REAL(GET_SLOT(x, Matrix_xSym));
169 :    
170 :     /* increment output count by number of off-diagonals */
171 :     nout = nnz;
172 :     for (j = 0; j < n; j++) {
173 :     int p2 = p[j+1];
174 :     for (jj = p[j]; jj < p2; jj++) {
175 :     if (iv[jj] != j) nout++;
176 :     }
177 :     }
178 :     SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
179 :     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nout));
180 :     ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
181 :     SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, nout));
182 :     aj = INTEGER(GET_SLOT(ans, Matrix_jSym));
183 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nout));
184 :     ax = REAL(GET_SLOT(ans, Matrix_xSym));
185 :     pos = 0;
186 :     for (j = 0; j < n; j++) {
187 :     int p2 = p[j+1];
188 :     for (jj = p[j]; jj < p2; jj++) {
189 :     int ii = iv[jj];
190 :     double xx = xv[jj];
191 : maechler 534
192 : bates 70 ai[pos] = ii; aj[pos] = j; ax[pos] = xx; pos++;
193 :     if (ii != j) {
194 :     aj[pos] = ii; ai[pos] = j; ax[pos] = xx; pos++;
195 :     }
196 :     }
197 :     }
198 :     UNPROTECT(1);
199 :     return ans;
200 :     }
201 : bates 153
202 : bates 478 SEXP dsCMatrix_ldl_symbolic(SEXP x, SEXP doPerm)
203 : bates 153 {
204 : bates 209 SEXP Ax, Dims = GET_SLOT(x, Matrix_DimSym),
205 :     ans = PROTECT(allocVector(VECSXP, 3)), tsc;
206 :     int i, n = INTEGER(Dims)[0], nz, nza,
207 :     *Ap, *Ai, *Lp, *Li, *Parent,
208 :     doperm = asLogical(doPerm),
209 :     *P = (int *) NULL, *Pinv = (int *) NULL;
210 : bates 153
211 : bates 209
212 : maechler 534 if (CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'L') {
213 : bates 209 x = PROTECT(ssc_transpose(x));
214 :     } else {
215 :     x = PROTECT(duplicate(x));
216 :     }
217 :     Ax = GET_SLOT(x, Matrix_xSym);
218 :     nza = length(Ax);
219 :     Ap = INTEGER(GET_SLOT(x, Matrix_pSym));
220 :     Ai = INTEGER(GET_SLOT(x, Matrix_iSym));
221 :     if (doperm) {
222 : bates 366 int *perm, *iperm = Calloc(n, int);
223 :    
224 : bates 209 SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, n));
225 :     perm = INTEGER(VECTOR_ELT(ans, 2));
226 : bates 366 ssc_metis_order(n, Ap, Ai, perm, iperm);
227 :     ssc_symbolic_permute(n, 1, iperm, Ap, Ai);
228 :     Free(iperm);
229 : bates 209 }
230 : bates 153 SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));
231 : bates 209 Parent = INTEGER(VECTOR_ELT(ans, 0));
232 : bates 478 SET_VECTOR_ELT(ans, 1, NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
233 : bates 209 tsc = VECTOR_ELT(ans, 1);
234 : maechler 534 SET_SLOT(tsc, Matrix_uploSym, mkString("L"));
235 :     SET_SLOT(tsc, Matrix_diagSym, mkString("U"));
236 : bates 209 SET_SLOT(tsc, Matrix_DimSym, Dims);
237 :     SET_SLOT(tsc, Matrix_pSym, allocVector(INTSXP, n + 1));
238 :     Lp = INTEGER(GET_SLOT(tsc, Matrix_pSym));
239 : bates 366 R_ldl_symbolic(n, Ap, Ai, Lp, Parent, P, Pinv);
240 : bates 209 nz = Lp[n];
241 :     SET_SLOT(tsc, Matrix_iSym, allocVector(INTSXP, nz));
242 :     Li = INTEGER(GET_SLOT(tsc, Matrix_iSym));
243 :     SET_SLOT(tsc, Matrix_xSym, allocVector(REALSXP, nz));
244 :     for (i = 0; i < nza; i++) REAL(Ax)[i] = 0.00001;
245 :     for (i = 0; i < n; i++) REAL(Ax)[Ap[i+1]-1] = 10000.;
246 : bates 366 i = R_ldl_numeric(n, Ap, Ai, REAL(Ax), Lp, Parent, Li,
247 : bates 209 REAL(GET_SLOT(tsc, Matrix_xSym)),
248 :     (double *) R_alloc(n, sizeof(double)), /* D */
249 : bates 366 P, Pinv);
250 : bates 209 UNPROTECT(2);
251 : bates 153 return ans;
252 :     }
253 : bates 184
254 : bates 478 SEXP dsCMatrix_metis_perm(SEXP x)
255 : bates 184 {
256 :     SEXP pSlot = GET_SLOT(x, Matrix_pSym),
257 :     ans = PROTECT(allocVector(VECSXP, 2));
258 :     int n = length(pSlot) - 1;
259 : maechler 534
260 : bates 184 SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));
261 : maechler 534 SET_VECTOR_ELT(ans, 1, allocVector(INTSXP, n));
262 : bates 184 ssc_metis_order(n,
263 :     INTEGER(pSlot),
264 :     INTEGER(GET_SLOT(x, Matrix_iSym)),
265 :     INTEGER(VECTOR_ELT(ans, 0)),
266 :     INTEGER(VECTOR_ELT(ans, 1)));
267 :     UNPROTECT(1);
268 :     return ans;
269 :     }
270 :    

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