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 296 - (view) (download) (as text)
Original Path: pkg/src/sscMatrix.c

1 : bates 10 #include "sscMatrix.h"
2 :    
3 :     SEXP sscMatrix_validate(SEXP obj)
4 :     {
5 :     SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
6 :     int *Dim = INTEGER(GET_SLOT(obj, Matrix_DimSym));
7 :     char *val;
8 :    
9 :     if (length(uplo) != 1)
10 :     return ScalarString(mkChar("uplo slot must have length 1"));
11 :     val = CHAR(STRING_ELT(uplo, 0));
12 :     if (strlen(val) != 1)
13 :     return ScalarString(mkChar("uplo[1] must have string length 1"));
14 :     if (toupper(*val) != 'U' && toupper(*val) != 'L')
15 :     return ScalarString(mkChar("uplo[1] must be \"U\" or \"L\""));
16 :     if (Dim[0] != Dim[1])
17 :     return ScalarString(mkChar("Symmetric matrix must be square"));
18 :     csc_check_column_sorting(obj);
19 :     return ScalarLogical(1);
20 :     }
21 :    
22 :     SEXP sscMatrix_chol(SEXP x, SEXP pivot)
23 :     {
24 : bates 257 SEXP pSlot = GET_SLOT(x, Matrix_pSym), xorig = x;
25 :     int *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
26 :     *Ap = INTEGER(pSlot),
27 :     *Lp, *Parent, info,
28 :     lo = toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'L',
29 :     n = length(pSlot)-1,
30 :     nnz, piv = asLogical(pivot);
31 : bates 10 SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("sscChol")));
32 : bates 257 int *Flag = Calloc(n, int), *Lnz = Calloc(n, int),
33 :     *P = (int *) NULL, *Pinv = (int *) NULL;
34 :     double *Ax;
35 : bates 10
36 : bates 257 if (lo) {
37 :     x = PROTECT(ssc_transpose(x));
38 :     Ai = INTEGER(GET_SLOT(x, Matrix_iSym));
39 :     Ap = INTEGER(GET_SLOT(x, Matrix_pSym));
40 :     }
41 :     SET_SLOT(val, Matrix_uploSym, ScalarString(mkChar("L")));
42 :     SET_SLOT(val, Matrix_diagSym, ScalarString(mkChar("N")));
43 :     SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
44 :     SET_SLOT(val, Matrix_ParentSym, allocVector(INTSXP, n));
45 :     Parent = INTEGER(GET_SLOT(val, Matrix_ParentSym));
46 :     SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, n + 1));
47 :     Lp = INTEGER(GET_SLOT(val, Matrix_pSym));
48 :     Ax = REAL(GET_SLOT(x, Matrix_xSym));
49 : bates 10 if (piv) {
50 : bates 257 SEXP trip = PROTECT(sscMatrix_to_triplet(x));
51 :     SEXP Ti = GET_SLOT(trip, Matrix_iSym);
52 : bates 10
53 : bates 257 /* determine the permutation with Metis */
54 :     Pinv = Calloc(n, int);
55 :     SET_SLOT(val, Matrix_permSym, allocVector(INTSXP, n));
56 :     P = INTEGER(GET_SLOT(val, Matrix_permSym));
57 :     ssc_metis_order(n, Ap, Ai, P, Pinv);
58 :     /* create a symmetrized form of x */
59 :     nnz = length(Ti);
60 :     Ai = Calloc(nnz, int);
61 :     Ax = Calloc(nnz, double);
62 :     Ap = Calloc(n + 1, int);
63 :     triplet_to_col(n, n, nnz, INTEGER(Ti),
64 :     INTEGER(GET_SLOT(trip, Matrix_jSym)),
65 :     REAL(GET_SLOT(trip, Matrix_xSym)),
66 :     Ap, Ai, Ax);
67 :     }
68 :     ldl_symbolic(n, Ap, Ai, Lp, Parent, Lnz, Flag, P, Pinv);
69 :     nnz = Lp[n];
70 : bates 10 SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));
71 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));
72 : bates 257 SET_SLOT(val, Matrix_DSym, allocVector(REALSXP, n));
73 :     info = ldl_numeric(n, Ap, Ai, Ax,
74 :     Lp, Parent, Lnz,
75 :     INTEGER(GET_SLOT(val, Matrix_iSym)),
76 :     REAL(GET_SLOT(val, Matrix_xSym)),
77 :     REAL(GET_SLOT(val, Matrix_DSym)),
78 :     (double *) R_alloc(n, sizeof(double)), /* Y */
79 :     (int *) R_alloc(n, sizeof(int)), /* Pattern */
80 :     Flag, P, Pinv);
81 :     if (info != n)
82 :     error("Leading minor of size %d (possibly after permutation) is indefinite",
83 :     info + 1);
84 :     Free(Flag); Free(Lnz);
85 :     if (piv) {
86 :     UNPROTECT(1);
87 :     Free(Pinv); Free(Ax); Free(Ai); Free(Ap);
88 :     }
89 :     UNPROTECT(lo ? 2 : 1);
90 :     return set_factorization(xorig, val, "Cholesky");
91 : bates 10 }
92 :    
93 :     SEXP sscMatrix_matrix_solve(SEXP a, SEXP b)
94 :     {
95 : bates 257 SEXP Chol = get_factorization(a, "Cholesky"), perm,
96 : bates 10 val = PROTECT(duplicate(b));
97 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
98 :     *bdims = INTEGER(getAttrib(b, R_DimSymbol)),
99 : bates 257 *Li, *Lp, j, n = adims[1], ncol = bdims[1], piv;
100 :     double *Lx, *D, *in = REAL(b), *out = REAL(val), *tmp = (double *) NULL;
101 : bates 10
102 :     if (!(isReal(b) && isMatrix(b)))
103 :     error("Argument b must be a numeric matrix");
104 :     if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
105 :     error("Dimensions of system to be solved are inconsistent");
106 : bates 296 if (Chol == R_NilValue) Chol = sscMatrix_chol(a, ScalarLogical(1));
107 : bates 257 perm = GET_SLOT(Chol, Matrix_permSym);
108 :     piv = length(perm);
109 :     if (piv) tmp = Calloc(n, double);
110 :     Li = INTEGER(GET_SLOT(Chol, Matrix_iSym));
111 :     Lp = INTEGER(GET_SLOT(Chol, Matrix_pSym));
112 :     Lx = REAL(GET_SLOT(Chol, Matrix_xSym));
113 :     D = REAL(GET_SLOT(Chol, Matrix_DSym));
114 :     for (j = 0; j < bdims[1]; j++, in += n, out += n) {
115 :     if (piv) ldl_perm(n, out, in, INTEGER(perm));
116 :     else Memcpy(out, in, n);
117 :     ldl_lsolve(n, out, Lp, Li, Lx);
118 :     ldl_dsolve(n, out, D);
119 :     ldl_ltsolve(n, out, Lp, Li, Lx);
120 :     if (piv) ldl_permt(n, out, Memcpy(tmp, out, n), INTEGER(perm));
121 : bates 10 }
122 : bates 257 if (piv) Free(tmp);
123 : bates 10 UNPROTECT(1);
124 :     return val;
125 :     }
126 :    
127 :     SEXP sscMatrix_inverse_factor(SEXP A)
128 :     {
129 : bates 257 return R_NilValue; /* FIXME: Write this function. */
130 : bates 10 }
131 :    
132 :     SEXP ssc_transpose(SEXP x)
133 :     {
134 :     SEXP
135 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sscMatrix"))),
136 :     islot = GET_SLOT(x, Matrix_iSym);
137 :     int nnz = length(islot),
138 :     *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),
139 :     *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
140 :    
141 :     adims[0] = xdims[1]; adims[1] = xdims[0];
142 :     if (toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'U')
143 :     SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));
144 :     SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, xdims[0] + 1));
145 :     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
146 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
147 :     csc_components_transpose(xdims[0], xdims[1], nnz,
148 :     INTEGER(GET_SLOT(x, Matrix_pSym)),
149 :     INTEGER(islot),
150 :     REAL(GET_SLOT(x, Matrix_xSym)),
151 :     INTEGER(GET_SLOT(ans, Matrix_pSym)),
152 :     INTEGER(GET_SLOT(ans, Matrix_iSym)),
153 :     REAL(GET_SLOT(ans, Matrix_xSym)));
154 :     UNPROTECT(1);
155 :     return ans;
156 :     }
157 : bates 70
158 :     SEXP sscMatrix_to_triplet(SEXP x)
159 :     {
160 :     SEXP
161 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tripletMatrix"))),
162 :     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 :    
192 :     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 209 SEXP sscMatrix_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 :     *Lnz = (int *) R_alloc(n, sizeof(int)),
210 :     *Flag = (int *) R_alloc(n, sizeof(int)),
211 :     *P = (int *) NULL, *Pinv = (int *) NULL;
212 : bates 153
213 : bates 209
214 :     if (toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'L') {
215 :     x = PROTECT(ssc_transpose(x));
216 :     } else {
217 :     x = PROTECT(duplicate(x));
218 :     }
219 :     Ax = GET_SLOT(x, Matrix_xSym);
220 :     nza = length(Ax);
221 :     Ap = INTEGER(GET_SLOT(x, Matrix_pSym));
222 :     Ai = INTEGER(GET_SLOT(x, Matrix_iSym));
223 :     if (doperm) {
224 :     int *perm;
225 :     SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, n));
226 :     perm = INTEGER(VECTOR_ELT(ans, 2));
227 :     ssc_metis_order(n, Ap, Ai, perm, Flag);
228 :     ssc_symbolic_permute(n, 1, Flag, Ap, Ai);
229 :     }
230 : bates 153 SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));
231 : bates 209 Parent = INTEGER(VECTOR_ELT(ans, 0));
232 :     SET_VECTOR_ELT(ans, 1, NEW_OBJECT(MAKE_CLASS("tscMatrix")));
233 :     tsc = VECTOR_ELT(ans, 1);
234 :     SET_SLOT(tsc, Matrix_uploSym, ScalarString(mkChar("L")));
235 :     SET_SLOT(tsc, Matrix_diagSym, ScalarString(mkChar("U")));
236 :     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 :     ldl_symbolic(n, Ap, Ai, Lp, Parent, Lnz, Flag, P, Pinv);
240 :     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 :     i = ldl_numeric(n, Ap, Ai, REAL(Ax), Lp, Parent, Lnz, Li,
247 :     REAL(GET_SLOT(tsc, Matrix_xSym)),
248 :     (double *) R_alloc(n, sizeof(double)), /* D */
249 :     (double *) R_alloc(n, sizeof(double)), /* Y */
250 :     (int *) R_alloc(n, sizeof(int)), /* Pattern */
251 :     Flag, P, Pinv);
252 :     UNPROTECT(2);
253 : bates 153 return ans;
254 :     }
255 : bates 184
256 :     SEXP sscMatrix_metis_perm(SEXP x)
257 :     {
258 :     SEXP pSlot = GET_SLOT(x, Matrix_pSym),
259 :     ans = PROTECT(allocVector(VECSXP, 2));
260 :     int n = length(pSlot) - 1;
261 :    
262 :     SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));
263 :     SET_VECTOR_ELT(ans, 1, allocVector(INTSXP, n));
264 :     ssc_metis_order(n,
265 :     INTEGER(pSlot),
266 :     INTEGER(GET_SLOT(x, Matrix_iSym)),
267 :     INTEGER(VECTOR_ELT(ans, 0)),
268 :     INTEGER(VECTOR_ELT(ans, 1)));
269 :     UNPROTECT(1);
270 :     return ans;
271 :     }
272 :    

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