Matrix r4655
Loading...
Searching...
No Matches
cholmod-etc.c
Go to the documentation of this file.
1#include "Mdefines.h"
2#include "idz.h"
3#include "cholmod-etc.h"
4
5cholmod_common c ;
6cholmod_common cl;
7
8cholmod_factor *M2CHF(SEXP obj, int values)
9{
10 cholmod_factor *L = (cholmod_factor *) R_alloc(1, sizeof(cholmod_factor));
11 memset(L, 0, sizeof(cholmod_factor));
12 SEXP dim = PROTECT(GET_SLOT(obj, Matrix_DimSym)),
13 type = PROTECT(GET_SLOT(obj, install("type"))),
14 perm = PROTECT(GET_SLOT(obj, Matrix_permSym)),
15 colcount = PROTECT(GET_SLOT(obj, install("colcount"))),
16 x = PROTECT(getAttrib(obj, Matrix_xSym));
17 L->n = INTEGER(dim)[0];
18 L->minor = L->n; /* FIXME: could be wrong for obj <- new(...) */
19 L->ordering = INTEGER(type)[0];
20 if (L->ordering != CHOLMOD_NATURAL)
21 L->Perm = INTEGER(perm);
22 else {
23 /* cholmod_check_factor allows L->Perm == NULL,
24 but cholmod_copy_factor does not test, so it segfaults ...
25 */
26 int n = (int) L->n, *Perm = (int *) R_alloc(L->n, sizeof(int));
27 for (int j = 0; j < n; ++j)
28 Perm[j] = j;
29 L->Perm = Perm;
30 }
31 L->ColCount = INTEGER(colcount);
32 L->is_super = INTEGER(type)[2];
33 if (L->is_super) {
34 L->is_ll = 1;
35 L->is_monotonic = 1;
36 SEXP super = PROTECT(GET_SLOT(obj, install("super"))),
37 pi = PROTECT(GET_SLOT(obj, install("pi"))),
38 px = PROTECT(GET_SLOT(obj, install("px"))),
39 s = PROTECT(GET_SLOT(obj, install("s")));
40 L->super = INTEGER(super);
41 L->pi = INTEGER(pi);
42 L->px = INTEGER(px);
43 L->s = INTEGER(s);
44 L->nsuper = LENGTH(super) - 1;
45 L->ssize = ((int *) L->pi)[L->nsuper];
46 L->xsize = ((int *) L->px)[L->nsuper];
47 L->maxcsize = INTEGER(type)[4];
48 L->maxesize = INTEGER(type)[5];
49 UNPROTECT(4);
50 } else {
51 L->is_ll = INTEGER(type)[1];
52 L->is_monotonic = INTEGER(type)[3];
53 if (values && x != R_NilValue) {
54 SEXP p = PROTECT(GET_SLOT(obj, Matrix_pSym)),
55 i = PROTECT(GET_SLOT(obj, Matrix_iSym)),
56 nz = PROTECT(GET_SLOT(obj, install("nz"))),
57 nxt = PROTECT(GET_SLOT(obj, install("nxt"))),
58 prv = PROTECT(GET_SLOT(obj, install("prv")));
59 L->p = INTEGER(p);
60 L->i = INTEGER(i);
61 L->nz = INTEGER(nz);
62 L->next = INTEGER(nxt);
63 L->prev = INTEGER(prv);
64 L->nzmax = ((int *) L->p)[L->n];
65 UNPROTECT(5);
66 }
67 }
68 L->itype = CHOLMOD_INT;
69 L->dtype = CHOLMOD_DOUBLE;
70 if (values && x != R_NilValue) {
71 switch (TYPEOF(x)) {
72 case CPLXSXP:
73 L->x = COMPLEX(x);
74 L->xtype = CHOLMOD_COMPLEX;
75 break;
76 case REALSXP:
77 L->x = REAL(x);
78 L->xtype = CHOLMOD_REAL;
79 break;
80 default:
81 ERROR_INVALID_TYPE(x, __func__);
82 break;
83 }
84 }
85 UNPROTECT(5);
86 return L;
87}
88
89cholmod_sparse *M2CHS(SEXP obj, int values)
90{
91 cholmod_sparse *A = (cholmod_sparse *) R_alloc(1, sizeof(cholmod_sparse));
92 memset(A, 0, sizeof(cholmod_sparse));
93 SEXP dim = PROTECT(GET_SLOT(obj, Matrix_DimSym)),
94 p = PROTECT(GET_SLOT(obj, Matrix_pSym)),
95 i = PROTECT(GET_SLOT(obj, Matrix_iSym)),
96 x = PROTECT(getAttrib(obj, Matrix_xSym));
97 A->nrow = INTEGER(dim)[0];
98 A->ncol = INTEGER(dim)[1];
99 A->p = INTEGER(p);
100 A->i = INTEGER(i);
101 A->nzmax = ((int *) A->p)[A->ncol];
102 A->stype = 0;
103 A->itype = CHOLMOD_INT;
104 A->xtype = CHOLMOD_PATTERN;
105 A->dtype = CHOLMOD_DOUBLE;
106 A->sorted = 1;
107 A->packed = 1;
108 if (values && x != R_NilValue) {
109 switch (TYPEOF(x)) {
110 case CPLXSXP:
111 A->x = COMPLEX(x);
112 A->xtype = CHOLMOD_COMPLEX;
113 break;
114 case REALSXP:
115 A->x = REAL(x);
116 A->xtype = CHOLMOD_REAL;
117 break;
118 default:
119 ERROR_INVALID_TYPE(x, __func__);
120 break;
121 }
122 }
123 UNPROTECT(4);
124 return A;
125}
126
127cholmod_dense *M2CHD(SEXP obj, int trans)
128{
129 cholmod_dense *A = (cholmod_dense *) R_alloc(1, sizeof(cholmod_dense));
130 memset(A, 0, sizeof(cholmod_dense));
131 SEXP dim = PROTECT(GET_SLOT(obj, Matrix_DimSym)),
132 x = PROTECT(GET_SLOT(obj, Matrix_xSym));
133 int m = INTEGER(dim)[0], n = INTEGER(dim)[1];
134 A->nrow = ((trans) ? n : m);
135 A->ncol = ((trans) ? m : n);
136 A->nzmax = A->nrow * A->ncol;
137 A->d = A->nrow;
138 A->dtype = CHOLMOD_DOUBLE;
139 switch (TYPEOF(x)) {
140 case CPLXSXP:
141 {
142 Rcomplex *px = COMPLEX(x);
143 if (!trans)
144 A->x = px;
145 else {
146 Rcomplex *py = R_Calloc(A->nzmax, Rcomplex);
147 ztranspose2(py, px, m, n);
148 A->x = py; /* NB: caller must do R_Free(A->x) */
149 }
150 A->xtype = CHOLMOD_COMPLEX;
151 break;
152 }
153 case REALSXP:
154 {
155 double *px = REAL(x);
156 if (!trans)
157 A->x = px;
158 else {
159 double *py = R_Calloc(A->nzmax, double);
160 dtranspose2(py, px, m, n);
161 A->x = py; /* NB: caller must do R_Free(A->x) */
162 }
163 A->xtype = CHOLMOD_REAL;
164 break;
165 }
166 default:
167 ERROR_INVALID_TYPE(x, __func__);
168 break;
169 }
170 UNPROTECT(2);
171 return A;
172}
173
174SEXP CHF2M(cholmod_factor *L, int values)
175{
176 if (L->itype != CHOLMOD_INT)
177 error(_("wrong '%s'"), "itype");
178 if (values && L->xtype != CHOLMOD_REAL && L->xtype != CHOLMOD_COMPLEX)
179 error(_("wrong '%s'"), "xtype");
180 if (values && L->dtype != CHOLMOD_DOUBLE)
181 error(_("wrong '%s'"), "dtype");
182 if (L->n > INT_MAX)
183 error(_("dimensions cannot exceed %s"), "2^31-1");
184 if (L->super) {
185 if (L->maxcsize > INT_MAX)
186 error(_("'%s' would overflow type \"%s\""),
187 "maxcsize", "integer");
188 } else {
189 if (L->n == INT_MAX)
190 error(_("n+1 would overflow type \"%s\""),
191 "integer");
192 }
193 if (L->minor < L->n) {
194 if (L->is_ll)
195 error(_("leading principal minor of order %d is not positive"),
196 (int) L->minor + 1);
197 else
198 error(_("leading principal minor of order %d is zero"),
199 (int) L->minor + 1);
200 }
201 char cl[] = ".CHM.....";
202 cl[0] = (!values) ? 'n' : ((L->xtype == CHOLMOD_COMPLEX) ? 'z' : 'd');
203 memcpy(cl + 4, (L->is_super) ? "super" : "simpl", 5);
204 SEXP obj = PROTECT(newObject(cl)),
205 dim = PROTECT(GET_SLOT(obj, Matrix_DimSym));
206 INTEGER(dim)[0] = INTEGER(dim)[1] = (int) L->n;
207 if (L->ordering != CHOLMOD_NATURAL) {
208 SEXP perm = PROTECT(allocVector(INTSXP, L->n));
209 Matrix_memcpy(INTEGER(perm), L->Perm, L->n, sizeof(int));
210 SET_SLOT(obj, Matrix_permSym, perm);
211 UNPROTECT(1);
212 }
213 SEXP type = PROTECT(allocVector(INTSXP, 6)),
214 colcount = PROTECT(allocVector(INTSXP, L->n));
215 INTEGER(type)[0] = L->ordering;
216 INTEGER(type)[1] = (L->is_super) ? 1 : L->is_ll;
217 INTEGER(type)[2] = (L->is_super) ? 1 : 0;
218 INTEGER(type)[3] = (L->is_super) ? 1 : L->is_monotonic;
219 INTEGER(type)[4] = (L->is_super) ? (int) L->maxcsize : 0;
220 INTEGER(type)[5] = (L->is_super) ? (int) L->maxesize : 0;
221 Matrix_memcpy(INTEGER(colcount), L->ColCount, L->n, sizeof(int));
222 SET_SLOT(obj, install("type"), type);
223 SET_SLOT(obj, install("colcount"), colcount);
224 if (L->is_super) {
225 SEXP super = PROTECT(allocVector(INTSXP, L->nsuper + 1)),
226 pi = PROTECT(allocVector(INTSXP, L->nsuper + 1)),
227 px = PROTECT(allocVector(INTSXP, L->nsuper + 1)),
228 s = PROTECT(allocVector(INTSXP, L->ssize));
229 Matrix_memcpy(INTEGER(super), L->super, L->nsuper + 1, sizeof(int));
230 Matrix_memcpy(INTEGER(pi), L->pi, L->nsuper + 1, sizeof(int));
231 Matrix_memcpy(INTEGER(px), L->px, L->nsuper + 1, sizeof(int));
232 Matrix_memcpy(INTEGER(s), L->s, L->ssize, sizeof(int));
233 SET_SLOT(obj, install("super"), super);
234 SET_SLOT(obj, install("pi"), pi);
235 SET_SLOT(obj, install("px"), px);
236 SET_SLOT(obj, install("s"), s);
237 UNPROTECT(4);
238 } else if (values) {
239 SEXP p = PROTECT(allocVector(INTSXP, L->n + 1)),
240 i = PROTECT(allocVector(INTSXP, L->nzmax)),
241 nz = PROTECT(allocVector(INTSXP, L->n)),
242 nxt = PROTECT(allocVector(INTSXP, L->n + 2)),
243 prv = PROTECT(allocVector(INTSXP, L->n + 2));
244 Matrix_memcpy(INTEGER(p), L->p, L->n + 1, sizeof(int));
245 Matrix_memcpy(INTEGER(i), L->i, L->nzmax, sizeof(int));
246 Matrix_memcpy(INTEGER(nz), L->nz, L->n, sizeof(int));
247 Matrix_memcpy(INTEGER(nxt), L->next, L->n + 2, sizeof(int));
248 Matrix_memcpy(INTEGER(prv), L->prev, L->n + 2, sizeof(int));
249 SET_SLOT(obj, Matrix_pSym, p);
250 SET_SLOT(obj, Matrix_iSym, i);
251 SET_SLOT(obj, install("nz"), nz);
252 SET_SLOT(obj, install("nxt"), nxt);
253 SET_SLOT(obj, install("prv"), prv);
254 UNPROTECT(5);
255 }
256 if (values) {
257 SEXP x;
258 R_xlen_t nx = (R_xlen_t) ((L->is_super) ? L->xsize : L->nzmax);
259 if (L->xtype == CHOLMOD_COMPLEX) {
260 PROTECT(x = allocVector(CPLXSXP, nx));
261 Matrix_memcpy(COMPLEX(x), L->x, nx, sizeof(Rcomplex));
262 } else {
263 PROTECT(x = allocVector(REALSXP, nx));
264 Matrix_memcpy(REAL(x), L->x, nx, sizeof(double));
265 }
266 SET_SLOT(obj, Matrix_xSym, x);
267 UNPROTECT(1);
268 }
269 UNPROTECT(4);
270 return obj;
271}
272
273SEXP CHS2M(cholmod_sparse *A, int values, char shape)
274{
275 cholmod_sparse *A_ = A;
276 if (A->itype != CHOLMOD_INT)
277 error(_("wrong '%s'"), "itype");
278 if (values && A->xtype != CHOLMOD_REAL && A->xtype != CHOLMOD_COMPLEX)
279 error(_("wrong '%s'"), "xtype");
280 if (values && A->dtype != CHOLMOD_DOUBLE)
281 error(_("wrong '%s'"), "dtype");
282 if (A->nrow > INT_MAX || A->ncol > INT_MAX)
283 error(_("dimensions cannot exceed %s"), "2^31-1");
284 if (!A->sorted)
285 cholmod_sort(A, &c);
286 if (!A->packed || A->stype != 0)
287 A = cholmod_copy(A, A->stype, 1, &c);
288 char cl[] = "..CMatrix";
289 cl[0] = (!values) ? 'n' : ((A->xtype == CHOLMOD_COMPLEX) ? 'z' : 'd');
290 cl[1] = shape;
291 int m = (int) A->nrow, n = (int) A->ncol, nnz = ((int *) A->p)[A->ncol];
292 R_xlen_t n1a = (R_xlen_t) n + 1;
293 SEXP obj = PROTECT(newObject(cl)),
294 dim = PROTECT(GET_SLOT(obj, Matrix_DimSym)),
295 p = PROTECT(allocVector(INTSXP, n1a)),
296 i = PROTECT(allocVector(INTSXP, nnz));
297 INTEGER(dim)[0] = m;
298 INTEGER(dim)[1] = n;
299 Matrix_memcpy(INTEGER(p), A->p, n1a, sizeof(int));
300 Matrix_memcpy(INTEGER(i), A->i, nnz, sizeof(int));
301 SET_SLOT(obj, Matrix_pSym, p);
302 SET_SLOT(obj, Matrix_iSym, i);
303 if (values) {
304 SEXP x;
305 if (A->xtype == CHOLMOD_COMPLEX) {
306 PROTECT(x = allocVector(CPLXSXP, nnz));
307 Matrix_memcpy(COMPLEX(x), A->x, nnz, sizeof(Rcomplex));
308 } else {
309 PROTECT(x = allocVector(REALSXP, nnz));
310 Matrix_memcpy(REAL(x), A->x, nnz, sizeof(double));
311 }
312 SET_SLOT(obj, Matrix_xSym, x);
313 UNPROTECT(1);
314 }
315 if (A != A_)
316 cholmod_free_sparse(&A, &c);
317 UNPROTECT(4);
318 return obj;
319}
320
321SEXP CHD2M(cholmod_dense *A, int trans, char shape)
322{
323 if (A->xtype != CHOLMOD_REAL && A->xtype != CHOLMOD_COMPLEX)
324 error(_("wrong '%s'"), "xtype");
325 if (A->dtype != CHOLMOD_DOUBLE)
326 error(_("wrong '%s'"), "dtype");
327 if (A->d != A->nrow) /* MJ: currently no need to support this case */
328 error(_("leading dimension not equal to number of rows"));
329 if (A->nrow > INT_MAX || A->ncol > INT_MAX)
330 error(_("dimensions cannot exceed %s"), "2^31-1");
331 int m = (int) A->nrow, n = (int) A->ncol;
332 if ((Matrix_int_fast64_t) m * n > R_XLEN_T_MAX)
333 error(_("attempt to allocate vector of length exceeding %s"),
334 "R_XLEN_T_MAX");
335 char cl[] = "...Matrix";
336 cl[0] = (A->xtype == CHOLMOD_COMPLEX) ? 'z' : 'd';
337 cl[1] = shape;
338 cl[2] = (shape == 'g')
339 ? 'e' : ((shape == 's') ? 'y' : ((shape == 'p') ? 'o' : 'r'));
340 SEXP obj = PROTECT(newObject(cl)),
341 dim = PROTECT(GET_SLOT(obj, Matrix_DimSym));
342 INTEGER(dim)[0] = (trans) ? n : m;
343 INTEGER(dim)[1] = (trans) ? m : n;
344 SEXP x;
345 if (A->xtype == CHOLMOD_COMPLEX) {
346 PROTECT(x = allocVector(CPLXSXP, (R_xlen_t) m * n));
347 Rcomplex *px = COMPLEX(x), *py = (Rcomplex *) A->x;
348 if (!trans)
349 Matrix_memcpy(px, py, (R_xlen_t) m * n, sizeof(Rcomplex));
350 else
351 ztranspose2(px, py, m, n);
352 } else {
353 PROTECT(x = allocVector(REALSXP, (R_xlen_t) m * n));
354 double *px = REAL(x), *py = (double *) A->x;
355 if (!trans)
356 Matrix_memcpy(px, py, (R_xlen_t) m * n, sizeof(double));
357 else
358 dtranspose2(px, py, m, n);
359 }
360 SET_SLOT(obj, Matrix_xSym, x);
361 UNPROTECT(3);
362 return obj;
363}
Matrix_int_fast64_t
long long Matrix_int_fast64_t
Definition Mdefines.h:27
_
#define _(String)
Definition Mdefines.h:44
ERROR_INVALID_TYPE
#define ERROR_INVALID_TYPE(_X_, _FUNC_)
Definition Mdefines.h:204
SET_SLOT
#define SET_SLOT(x, what, value)
Definition Mdefines.h:86
GET_SLOT
#define GET_SLOT(x, what)
Definition Mdefines.h:85
newObject
SEXP newObject(const char *)
Definition objects.c:4
Matrix_permSym
SEXP Matrix_permSym
Definition Msymbols.h:18
Matrix_DimSym
SEXP Matrix_DimSym
Definition Msymbols.h:3
Matrix_xSym
SEXP Matrix_xSym
Definition Msymbols.h:22
Matrix_iSym
SEXP Matrix_iSym
Definition Msymbols.h:13
Matrix_pSym
SEXP Matrix_pSym
Definition Msymbols.h:17
CHD2M
SEXP CHD2M(cholmod_dense *A, int trans, char shape)
Definition cholmod-etc.c:321
M2CHF
cholmod_factor * M2CHF(SEXP obj, int values)
Definition cholmod-etc.c:8
M2CHS
cholmod_sparse * M2CHS(SEXP obj, int values)
Definition cholmod-etc.c:89
CHF2M
SEXP CHF2M(cholmod_factor *L, int values)
Definition cholmod-etc.c:174
c
cholmod_common c
Definition cholmod-etc.c:5
cl
cholmod_common cl
Definition cholmod-etc.c:6
CHS2M
SEXP CHS2M(cholmod_sparse *A, int values, char shape)
Definition cholmod-etc.c:273
M2CHD
cholmod_dense * M2CHD(SEXP obj, int trans)
Definition cholmod-etc.c:127
trans
void trans(SEXP p0, SEXP i0, SEXP x0, SEXP p1, SEXP i1, SEXP x1, int m, int n)
Definition coerce.c:3356
Matrix_memcpy
void * Matrix_memcpy(void *dest, const void *src, R_xlen_t length, size_t size)
Definition utils.c:70
Generated by doxygen 1.9.8